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FITTING  SURFACES  TO  SCATTERED  DATA 


Larry  L.  Schumaker 


This  paper  Is  a  survey  of  a  variety  of  numerical  methods 
for  fitting  a  function  to  data  given  at  a  set  of  points  scat¬ 
tered  throughout  a  domain  in  the  plane.  We  discuss  four 
classes  of  methods:  (1)  global  interpolation,  (2)  local  inter¬ 
polation,  (3)  global  approximation,  and  (4)  local  approximation. 
W«c.also  discus^' two-stage  methods  and  contouring.  The  surfaces 
constructed  will  include  polynomials,  spline  functions,  and  ra¬ 
tional  functions,  among  others. 

1.  Introduction 

Our  aim  is  to  survey  methods  for  solving  the  following 
problem. 

PROBLEM  1.1.  Let  D  be  a  domain  in  the  (x.y)- plane,  and  sup¬ 
pose  F  is  a  real-valued  function  defined  on  D.  Suppose  we 
are  given  the  values  Fj^  =  F(x^, y^)  of  F  at  some  set  of 
points  (Xj^,y^)  located  in  D,  i  =  1, 2,  ...,N.  Find  a  function 
f 


y^)  located  in 
defined  on  D  which  reasonably  approximates  F. 


This  problem  is,  of  course,  precisely  the  problem  of  fit¬ 
ting  a  surface  to  given  data.  In  many  cases  the  domain  D  is 
a  rectangle  and  the  data  points  lie  on  a  rectangular  grid. 

There  are,  however,  many  practical  problems  (see  the  following 
section  for  some  specific  examples),  where  D  is  of  unusual 
shape  and  where  the  data  points  are  irregularly  scattered 
throughout  D.  Thus,  while  we  shall  pay  some  attention  to  spe¬ 
cial  methods  for  regularly  spaced  data,  we  are  actually  more 
interested  in  the  general  case. 

There  are  basically  two  approaches  to  handling  Problem  1.1. 

First,  we  may  try  to  construct  a  function  f  which  interpolates 
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the  data  exactly;  i.e.,  such  that 

(1.1)  f(x^,y^)  =  i  =  1,2,..,,N. 

This  approach  may  be  desirable  when  the  function  values  at  the 
data  points  are  known  to  high  precision  and  where  it  is  highly 
desirable  that  these  values  be  preserved  by  the  approximating 
function. 

The  second  approach  involves  constructing  f  which  only 
approximately  fits  the  data.  This  may  be  regarded  as  data 
smoothing  and  will  be  desirable  when  (as  is  often  the  case) 
the  data  are  subject  to  inaccurate  measurement  or  even  errors. 
The  question  of  whether  interpolation  or  approximation  should 
be  used  will  not  be  discussed  further  here--this  is  a  problem 
which  must  be  settled  for  the  individual  problem  at  hand. 

In  discussing  Problem  1.1,  it  will  be  convenient  to  make 
a  further  distinction  between  those  methods  which  are  local  in 
character  (i.e.,  where  the  value  of  the  constructed  surface  f 
at  the  point  (x, y)  depends  only  on  the  data  at  relatively 
nearby  points)  and  those  methods  which  are  global  in  nature. 
Thus,  we  discuss  four  categories  of  methods  in  sections  3-6: 

(1)  global  interpolation,  (2)  local  interpolation,  (3)  global 
approximation,  and  (4)  local  approximation.  In  each  of  these 
sections  we  further  subdivide  the  material  according  to  the 
type  of  functions  being  used  and  the  type  of  data  (scattered  or 
not)  for  which  the  method  is  suitable. 

In  discussing  methods  which  apply  only  to  special  arrange¬ 
ments  of  data  points,  we  have  two  objectives  in  mind.  First, 
the  methods  are  of  interest  in  their  own  right.  More  important¬ 
ly  in  terms  of  Problem  1.1,  however,  such  methods  can  also  be 
used  in  two-stage  processes  in  which  we  first  construct  a  sur¬ 
face  g  based  on  the  scattered  data,  and  then  use  g  to  gen¬ 
erate  regular  data  for  the  construction  of  another  (perhaps 
smoother  or  more  convenient)  surface  f.  Such  two- stage  methods 
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will  be  discussed  (along  with  several  examples)  in  more  detail 
in  section  7. 

For  many  of  the  methods  based  on  regular  data  and  some  of 
those  for  scattered  data,  error  bounds  are  available  to  indi¬ 
cate  how  well  smooth  functions  are  approximated  by  the  surface 
constructed.  We  do  not  have  space  to  go  into  the  extensive 
literature  on  error  bounds,  A  simple  test  of  how  well  a  method 
will  approximate  smooth  functions  is,  however,  provided  by  its 
ability  to  reproduce  polynomial  surfaces  exactly  (that  is,  if 
F  is  a  polynomial  in  x  and  y  up  to  a  certain  degree,  then 
the  surface  f  is  identically  equal  to  F) .  For  many  of  the 
methods  we  will  be  able  to  indicate  the  corresponding  degree  of 
exactness. 

In  many  of  the  applications  of  surface-fitting  techniques 
(cf.  the  examples  in  section  2),  the  ultimate  aim  is  to  use  the 
data  to  construct  a  contour  map  of  the  unknown  function.  Since 
F  is  known  only  at  the  data  points,  we  must  be  '~ontent  to  con¬ 
struct  a  contour  map  for  one  of  our  fitted  surfaces.  In  sec¬ 
tion  8  we  discuss  some  approaches  to  accomplishing  this  numeri¬ 
cally. 

We  close  this  introduction  with  a  disclaimer-- this  survey 
does  not  include  all  possible  methods  for  fitting  surfaces  to 
scattered  data.  For  example,  we  have  not  discussed  Fourier 
series  methods,  spatial  filtering,  and  other  such  related  sta¬ 
tistical  techniques.  In  addition,  the  set  of  references  for 
those  methods  which  we  have  discussed  are  also  not  complete. 

My  original  intention  was  to  compile  as  complete  a  bibliography 
as  possible,  but  the  sheer  bulk  of  relevant  papers  and  my  in¬ 
ability  to  locate  all  of  them  convinced  me  to  settle  for  less. 

I  have  opted  to  quote  a  fairly  representative  list  of  papers, 
including  several  other  surveys.  Further  references  can  be 
found  by  consulting  these.  I  shall  be  very  happy  to  receive 
information  on  references  and  methods  I  have  overlooked. 


A 
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2 .  Examples 


In  this  section  we  shall  quote  several  explicit  examples 
of  Problem  1.1  to  emphasize  the  fact  that  unusually  shaped  re¬ 
gions  and  scattered  data  do  arise  frequently  in  practice, 

EXAMPLE  2.1.  Petroleum  exploration.  In  exploring  for  petro¬ 
leum,  the  contours  of  various  underground  layers  of  sandstone, 
shale,  limestone,  etc.  can  be  important  indicators  of  possible 
oil  fields.  Frequently,  data  on  such  layers  is  available  from 
exploratory  wells,  which,  however,  have  most  likely  been  drilled 
at  locations  scattered  randomly  throughout  some  geographical  re¬ 
gion  of  interest.  To  quote  a  specific  example,  Robinson, 
Charlesworth,  and  Ellis  [166]  consider  precisely  this  problem 
for  some  data  obtained  from  7,500  wells  drilled  in  Alberta.  For 
another  example  of  this  type,  see  Whitten  and  Koelling  [208]. 

Problems  similar  to  that  mentioned  in  Example  2.1  arise 
frequently  in  cartography  and  submarine  topography  where  the 
measurements  represent  actual  elevations.  In  some  cases  the 
measurements  must  be  taken  from  photographs  or  from  sonar  mea¬ 
surements  and  are  usually  subject  to  some  measurement  error 
(eg.  see  Kubik  [125]  for  a  discussion  of  photogrammetry)  . 

EXAMPLE  2.2.  Geological  maps.  There  are  a  great  many  problems 
in  Geology  and  the  earth  sciences  in  which  the  data  arises  from 
some  other  function  of  location  besides  actual  elevations.  For 
example,  some  geological  variables  of  interest  might  include 
concentrations  of  various  chemicals,  specific  gravity,  electri¬ 
cal  resistivity,  grain  size,  texture,  optical  properties,  iso¬ 
tope  ratios,  etc.  To  quote  a  specific  example,  Bhattacharyya 
[21,  22]  discusses  methods  for  fitting  a  surface  to  measurements 
(taken  by  airborne  sensors)  of  magnetic  potentials  over  a  cer¬ 
tain  portion  of  the  Yukon.  See  also  Bhattacharyya  and  Raychaud- 
huri  [23]  and  Crain  and  Bhattacharyya  [61]. 
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The  importance  of  surface-fitting  methods  in  the  earth  sci¬ 
ences  can  be  judged  by  the  large  number  of  papers  in  the  area 
relating  to  various  fitting  methods.  For  a  further  list  of 
problems  and  a  discussion  of  some  of  the  methods  which  have  been 
applied,  see  the  books  of  Bohrenberg  and  Giese  [31],  Chorley 
[51],  David  [62],  Harbaugh  and  Merriam  [98],  and  Merriam  [140], 
Recent  survey  papers  induce  Whitten  [20  3  ,  205]  and  Whittfi.-i  and 
Koelling  [20  7],  To  add  just  a  few  more  of  the  papers  in  the 
geological  literature  dealing  with  surface  fitting  tc  our  list, 
we  mention  Anderson  [7],  Grant  [91],  Messing,  Lee,  and  Pierce 
[114],  Holroyd  and  Bhattacharyya  [115],  Kubik  [123,  125],  Nor- 
cliffe  [151],  Reilly  [162],  Whitten  [200,  201,  204],  and  Whit¬ 
ten  and  Koelling  [206], 


EXAMPLE  2.3.  Heart  potentials.  In  order  to  diagnose  certain 
abnormal  heart  conditions,  it  is  desired  to  make  a  series  of 
several  hundred  contour  maps  of  the  heart  potential  field  at 
time  steps  of  1/100  of  a  second  throughout  a  heart  beat.  Data 
on  these  heart  potr ntials  can  be  obtained  by  fitting  the  patient 
with  a  shirt  containing  probes.  Because  of  body  geometry,  when 
this  shirt  is  flattened  out  it  takes  the  nonrectangular  form 
illustrated  in  Figure  1,  Although  the  probes  could  be  arranged 
fairly  regularly  in  this  domain,  because  of  the  added  signifi- 


*  *  •  .  «  • 
.  •  • 

*  •  I  •  •  • 


Figure  1,  Heart  Potential  Measurements 
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cance  of  frontal  measurements,  in  practice  more  probes  are 
fitted  there  than  in  the  back.  This  example  was  brought  to  my 
attention  by  Ms.  Patrizia  Ciarlini  of  Rome. 

Potential  fields  arise  in  many  other  applications.  We 
have  already  mentioned  Geology  in  Example  2.2.  For  some  exam¬ 
ples  in  modelling  plasmas  see  Buneman  [40].  The  problem  arises 
in  Biersack  and  Fink  [24]  in  experimentally  studying  crystal 
structure  using  neutron  bombardment.  Data  from  waveform  dis¬ 
tortion  in  electronic  circuits  can  be  found  in  Akima  [5,  6]. 


3 .  Global  interpolation  methods 

In  this  section  we  outline  several  methods  for  solving  the 
interpolation  problem  (1.1). 

3.1  Polynomial  interpolation.  (Scattered  data) .  The  general 
theory  of  finite  dimensional  interpolation  is,  of  course,  very 
well  known  (e.g.,  see  Davis  [63]).  Briefly,  if  N 

functions  defined  on  the  domain  D,  then  the  function 


N 

(3.1)  f(x,y)  -  Ea.0.(x,y) 

j=l  J  J 


will  satisfy  (1.1)  if  and  only  if 
linear  system 
n 


(3.2) 


j  =  l 


a.0.(x^, 


i 


(a.i“ 


is  a  solution  of  the 


1..2, 


•  f 


N. 


This  system  has  a  (unique)  solution  for  arbitrary  choices  of 
data  precisely  when  it  is  nonsingular.  This  depends  on  the 

choice  of  functions  {0.}^  and  the  location  of  the  data  points. 

J 

To  illustrate  this  method,  we  may  choose  the 

polynomials  in  x  and  y.  Given  N,  there  is  some  leeway  in 

the  choice  of  which  powers  of  x  and  y  to  use.  For  example, 

with  N  =  3  one  could  use  the  functions  1,  x,  y  or  possibly 

2  2 

the  functions  1,  x  ,  y  ,  etc.  When  N  is  of  the  form  N  = 


■j 


(d  +  l)(d-tl),  we  might  use  the  functions 


y) 


r  V  u, 

lx  y  J 


d 

VvO 


0 


As  simple  as  this  sounds,  there  are  some  serious  difficul¬ 
ties  with  polynomial  interpolation  of  scattered  data.  For  open¬ 
ers,  it  is  not  so  easy  to  guarantee  that  the  system  (3.2)  is 
nonsingular.  To  give  a  very  simple  example,  consider  the  case 
N  =  3  with  the  functions  1,  x,  y.  If  the  three  data  points 
happen  to  lie  on  a  line,  then  (3.2)  will  in  fact  be  singular. 
Even  when  (3.2)  is  nonsingular,  it  will  often  be  the  case  (at 
least  if  N  is  moderately  large)  that  the  system  will  be  ill- 
conditioned.  Finally,  as  is  well  known,  polynomials  of  even 
moderate  degree  exhibit  a  considerable  oscillatory  character, 
and  the  resulting  surface  (even  though  it  is  is  often  too 
undulating  to  be  acceptable.  The  general  problem  of  polynomial 
interpolation  to  scattered  data  is  not  usually  treated  in  Nu¬ 
merical  Analysis  and  Approximation  Theory  books  (see,  however, 
Kunz  [126],  Prenter  [157],  and  Steffenson  [186]).  Some  papers 
dealing  with  the  question  include  Guenther  [93],  Thatcher  [189], 
Thatcher  and  Milne  [190],  and  Whaples  [197],  Assuming  the  in- 
terpolant  exists,  error  bounds  have  been  studied  in  Ciarlet  and 
Raviart  [52-55]. 

Let 


(3.3) 


P 

m,  n 


span 


lx 


V  fi-im  ,n 

y  ^=0,^=0 


be  the  space  of  polynomials  of  degree  m  in  x  and  of  degree 
n  in  y.  This  linear  space  is  of  dimension  (m+1) (n+1)  and 
is,  in  fact,  the  tensor  product  of  the  linear  spaces  and 

P^.  It  is  perhaps  of  interest  to  note  that  there  always  exists 
a  (usually  nonunique)  polynomial  p  e  ^  which  solves  the 

IN 

interpolation  problem  (1.1),  no  matter  how  the  data  points  are 
positioned,  see  Prenter  [158]. 


3,2  Polynomial  interpolation  (^ridded  data) .  We  begin  this 
subsection  by  defining  what  we  mean  by  gridded  data.  Let 

(3.4)  H  -  [a, b]  X  [c,d] 


be  a  rectangle,  and  let 


(3.5) 


a 

c 


<  X 


k+1 


b 

d. 


We  suppose  now  that  F 
we  have  the  values  of 
lar  grid  defined  by  (3, 

(3.6)  F^,  F(x.,y^), 


is  a  function  defined  on  H,  and  that 
F  at  the  corner  points  of  the  rectangu- 
5);  i.e., 

i  =  0,1,..  ,,k+l 
j  0,1,..., i+l. 


This  is  a  total  of  N  =  (k+2)  (.£+2)  data  points. 

It  is  quite  easy  to  show  that  there  exists  a  unique  poly 
nomial  p  in  the  class  ^^-,01  definition  (3.3)) 

which  interpolates  the  gridde<l  data  given  in  (3.4)- (3.6).  In 
fact,  p  can  be  written  down  explicitly  in  terms  of  the  one¬ 
dimensional  Lagrange  polvnomials  as 


k+1  i+1 

(3,7)  p(x,y)  =  L  Z  F  L  (x)L  (y), 
i=0  j=0  J  ^ 

where  the  {L^(x)]q^^  and  (Lj(y)]Q^^  are  the  usual  one¬ 
dimensional  Lagrange  polynomials  associated  with  the  interpola¬ 
tion  points  respectively.  Interpolation 

of  gridded  data  by  polynomials  has  been  discussed  in  various 
books  and  papers--we  do  not  bother  with  a  long  list.  See  e.g, 
Prenter  [157]  or  Steffenson  [186],  More  recently,  there  has 
been  considerable  work  on  Hermite  and  osculatory  interpolation 
in  several  variables;  see  e.g.  Ahlin  [3],  Haussman  [99,101,102], 
and  Salzer  [168-170]. 
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3,3  Shepard's  method.  In  this  subsection  we  discuss  a  method 
of  Shepard  [180]  and  some  modifications  of  it.  The  method  ap¬ 
plies  to  arbitrarily  spaced  data,  and  the  interpolating  function 
can  be  written  down  explicitly. 

Let  p  be  some  metric  in  the  plane,  for  example  the  usual 
distance  metric.  Given  a  point  (x, y),  let  =  p( (x, y) , (x^, y^)) 


for  i  =  1,  2,  .  .  ,  ,N.  Let  0  <  p  <  oo. 
tion  formula  is  defined  by 

r 


(3.8)  f(x,y)  =  { 


Then  Shepard's  interpola- 


when  r^  ^  0,  all  i 


when  r . 


0. 


The  formula  (3.8)  is  defined  for  all  points  (x, y)  in  the 
2 

plane  R  .  It  is  clear  ..rom  the  definition  that  it  interpolates 
the  values  at  the  data  points  (x^, y^),  i  =  1,2, ...,N. 

The  value  of  f(x, y)  at  nondata  points  is  obtained  as  a  weight¬ 
ed  average  of  all  the  data  values,  where  the  i*"^  measurement 
is  weighted  according  to  the  distance  of  (x,y)  from  the  point 
(Xj_,y^)  . 

We  shall  briefly  recount  some  of  the  properties  of  Shep¬ 
ard's  formula.  First,  by  converting  all  of  the  terms  to  a 
common  denominator,  it  can  be  shown  that 
N 

(3.9)  f(x,y)  =  Z  F/.(x,y), 

i=l  ^  ^ 

where 

TT  Ir.Cx.y)!** 

j  =  l  J 

(3.10)  Aj(*,y)  - 

Z  Z  [r  (x,y)]^ 
k=l  ^ 

These  functions  satisfy 


i  =  1,  2, . .  .,N  . 
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(3.11)  A^(x^,yj)  -  “  1,2, 

The  representation  (3.9)  is  numerically  more  stable  than  the 
original  formula  (3.8). 

In  view  of  its  definition,  we  see  that  the  function  f(x, y) 
constructed  by  Shepard  is  not  a  simple  polynomial  or  rational 
function.  It  is  clear,  however,  that  except  for  the  points 
(x^fYi),  it  is  analytic  everywhere  in  the  plane.  Its  behavior 
in  the  vicinity  of  the  data  points  depends  on  the  size 

of  n.  It  can  be  shown  that  for  0  <  p  <  1,  f  has  cusps  at 
these  points.  For  1  <  ^,  f  has  flat  spots  at  the  data  points 
(i.e.,  the  partial  derivatives  vanish  there).  We  also  observe 
the  interesting  property  that 

(3.12)  min  F  <  f(x, y)  <  max  F  . 

ISisN  ^  liigN  ^ 

We  may  also  note  that  if  the  data  came  from  a  constant  function, 
i.e.,  F^  =  c,  i  =  1,2, ...,N,  then  f  is  also  the  constant 
function  f  ■  c. 

We  now  comment  on  the  choice  of  p.  To  get  smooth  surfaces 
without  cusps,  it  is  desirable  to  take  1  <  p.  On  the  other 
hand,  if  p  is  relatively  large,  then  the  surface  tends  to  be¬ 
come  very  flat  near  the  data  points  and  consequently  quite  steep 
at  points  in  between.  Experiments  (cf.  Gordon  and  Wixom  [90], 
Poeppelmeir  [155],  and  Shepard  [180])  seem  to  indicate  that  a 
choice  of  p  =  2  is  perhaps  a  good  tradeoff.  ([155]  contains 
several  examples  showing  the  behavior  as  a  function  of  p.) 

There  are  several  drawbacks  to  Shepard's  method  (3.8),  as 
pointed  out  by  Shepard  [180]  himself.  First,  if  N  is  large, 
then  there  is  a  very  considerable  amount  of  calculation  in¬ 
volved  in  evaluating  f(x,y)  at  a  particular  point.  Secondly, 
the  weights  are  assigned  on  the  basis  of  the  distance  of  points 
from  (x, y)  only,  not  their  direction.  Finally,  the  flat  spots 


4 


in  the  neighborhood  of  the  data  points  is  somewhat  disturbing. 
The  first  of  these  objections  can  be  met  by  defining  a  local 
version  of  the  formula,  which  we  shall  do  in  section  4,5.  It  is 
possible  to  construct  an  analogous  fonnala  which  accounts  for 
direction.  For  details,  see  Shepard  [180].  Finally,  we  briefly 
discuss  handling  the  flat  spots. 

Suppose  in  addition  to  the  function  values  F^  at  each 
point  (x^jy^)  we  also  have  estimates 


FX^  and  FY^ 


of 


N 


Then  we  may  consider  the  function 


(3.13)  f(x,y)  =  EA.(x,y)[F^  +  (x-x^)FX.  +  (y-y^)FYj, 

i=l 


It  is  easily  checked  that  this  function  also  interpolates,  and 
that 


(3.14)  =  FX^,  fy(vyi>  =  ^^’i'  i  =  1,2,...,N. 

This  property  may  be  expressed  in  the  assertion  that  if  the 
data  came  from  a  plane  surface  F,  then  f  will 

exactly  reproduce  this  surface.  To  use  formula  (3,13)  in  prac¬ 
tice  on  the  data- fitting  Problem  1.1,  we  have  to  carry  out  a 
two-stage  approximation  process  in  which  the  first  stage  con¬ 
sists  of  some  method  for  estimating  the  slope  at  each  of  the 
data  points. 

It  might  be  of  practical  interest  in  some  cases  to  con¬ 
struct  still  a  more  sophisticated  version  of  Shepard's  formula 
which  would  exactly  reproduce  higher-order  polynomial  surfaces. 
One  approach  to  doing  this  is  to  use  the  following  lemma. 

LEMMA  3.1.  (Barnhill  [15]).  Let  P  and  Q  be  linear  projec¬ 
tions  of  some  linear  space  of  functions  ^  into  itself.  Sup¬ 
pose  that  Q  exactly  reproduces  the  linear  subspace  E  c:  y- 
i.e. . 

(3.15)  Qp  =  p,  all  p  e  E. 
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In  addition,  suppose  that  is  a  set  of  linear  functionals 

on  and  that 

(3 . 16)  all  £  e  i  =  1,  2,  . .  .  ,m. 

Then  the  Boolean  sum  projector 

(3.17)  P  ©  Q  ■=  P  +  Q  -  PQ 

enjoys  the  function  precision  of  Q  (i.e.,  reproduces  E)  and 
the  interpolation  properties  of  P  (i.e.,  (3.16)  also  holds 
for  P  ©  Q)  . 

This  result  permits  the  construction  of  interpolation 
schemes  using  Shepard's  formula  which  reproduce  higher-order 
surfaces.  For  an  example,  see  Poeppelmeir  [155]  where  Shepard's 
formula  is  combined  with  a  certain  local  interpolation  scheme 
which  reproduces  quadratic  surfaces.  In  closing  this  section 
we  note  that  Shepard's  formula  can  also  be  interpreted  as  aris¬ 
ing  from  weighted  least  squares--see  section  5.1. 

3.4  Spline  interpolation  (scattered  data) .  Suppose  X  is  a 
linear  space  of  "smooth"  functions  defined  on  the  domain  D, 
and  let 

(3.18)  U  =  {f  e  X:  f(x.,y^)  =  F^,  i  =  1,2,  ...,N}. 

U  is  the  set  of  smooth  functions  which  interpolate.  Now  sup¬ 
pose  that  0  is  a  functional  on  X  which  measures  the  smooth¬ 
ness  of  an  element  in  X--the  smaller  0(f)  is,  the  smoother 
f  is.  Then  we  may  consider  the  following  minimization  problem: 

(3.19)  Find  s  e  U  such  that  0(s)  =  inf  0(u). 

ueU 

The  function  s  will  be  the  smoothest  interpolant,  and  in  view 
of  the  similarity  with  classical  spline  approximation,  s  is 
called  a  spline  function  Interpolating  F.  The  basic  questions 
concerning  spline  interpolation  center  around  existence,  unique¬ 
ness,  characterization,  and  construction.  A  quite  general 
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abstract  theory  of  spline  Interpolation  has  been  built  up  (see 
eg.  Laurent  [127]  and  references  therein).  In  this  section  we 
quote  some  specific  examples  which  can  be  used  on  Problem  1.1, 
Where  X  is  a  semi-Hilbert  space,  0(f)  =  ||f||,  where  H'll 
is  a  seminorm  on  X,  and  N  -  (f  e  X:  ||f||  =  O),  it  is  possible 
to  show  (under  some  additional  mild  conditions  on  X,  see  Duchon 
[72,  73])  that  problem  (3.19)  always  has  a  solution  which  is 
unique  up  to  an  element  in  N.  Moreover,  it  can  be  shown  that 
there  exists  a  reproducing  kernel  K  defined  on  DxD  such 
that 

N  d 

(3.20)  s(x,y)  =  E  a  K((x,y)  ;  (x  ,y  ))  +  E  b  p  (x,  y) , 

i=l  ^  ^  i=l 

,  -d 

where  is  a  basis  for  N.  Moreover,  the  coefficients 

[a^]  and  (b^)  can  be  determined  from  the  linear  system  of 
equations 

N  d 

(3.21)  E  K((x  ,y  )  ;  (x  ,y  ))a  +  Ebp  (x,y)=F,  j=l,...,N 

JJ  ^  i“l  ^  ^  J  J  J 

N 

E  a  p  (x.,y  )  -  0,  k  =  l,2,,.,,d, 

i^l  1  K  1  1 


The  development  with  semi-Hilbert  spaces  in  Duchon  [72,73] 
is  an  extension  of  earlier  work  of  Atteia  [10-12]  and  Thomann 
[192-193]  using  Hilbert  spaces.  The  essential  difficulty  in 
applying  the  general  result'  is  the  construction  of  an  appro¬ 
priate  reproducing  kernel.  We  turn  now  to  some  specific  exam¬ 
ples. 

Suppose  X  is  the  space  of  all  functions  on  the  rectangle 

D  =  H  (cf.  (3.4))  which  have  (distributional)  derivatives  up  to 

2 

order  2  which  lie  in  L  (H)  .  For  f  e  X,  let 


(3.22)  e(f)  =  ff  iD^f 

D  ^ 


+  2|D  D  fl^  +  iD^fl^. 

X  y  y 


2 
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The  reproducing  kernel  in  this  case  can  be  written  down  as  an 

infinite  series  involving  sin  and  cos,  and  the  space  N  is 

spanned  by  1,  x,  and  y.  Similarly,  if  we  replace  H  by  the 

unit  disc  UD,  the  kernel  can  be  computed  as  an  infinite  series 

(see  Atteia  [10-12]  and  Thomann  [192-193]).  Thomann  considers 

computation  of  these  splines  by  approximating  the  infinite 

series- -FORTRAN  programs  are  also  included. 

If  we  replace  the  bounded  sets  H  or  UD  by  the  entire 
2 

plane  R  and  introduce  an  appropriate  space  X,  it  is  possible 
to  obtain  explicit  expressions  for  the  reproducing  kernel.  This 

^3 

is  the  content  of  Duchon  [72, /3].  In  particular,  let  H 

2 

be  the  set  of  all  tempered  distributions  f  on  R  whose 
Fourier  transforms  f  satisfy  /jfjt^^dt  <  »,  Let  de¬ 

note  the  set  of  all  functions  which  have  derivatives  up  to 

<^3 

order  m  lying  in  H  ,  Our  first  example  concerns  the  space 
20 

X  .  If  we  choose  0  as  in  (3.22),  then  the  interpolating 
spline  solution  of  (3.19)  is  of  the  form 
N  2 

(3.23)  s(x,y)  =  Z  a^r^(x,y)  log  (r^(x,y))  +  b^^x  +  b^y  +  b^, 

where  r^(x, y)  =  [(x-x^)^+  (y-y^)  ]^.  The  coefficients  are  de¬ 
termined  from  the  system  (3.21)  with  d  =  3,  N  =  span  {l,x,y), 
and  K(z,w)  =  |z-w|  log(z-w|.  Duchon  refers  to  this  type  of 
spline  as  a  thin  plate  spline  since  the  expression  0  relates 

to  the  energy  in  a  thin  plate  forced  to  interpolate  the  data. 

2 

This  spline  belongs  to  C(R  ). 

21 

As  a  second  example,  suppose  we  consider  X  =  X  .  In 
this  case  the  solution  of  (3.19)  with  0  given  by  (3.22)  has 
the  form 

N 

(3.24)  s(x,y)  =  Z  a^(r^(x,y))  +  bj^x  +  b^y  +  b^. 

|z-w|  .  Duchon  [72,73]  refers  to  these  splines 


Here  K(z,w)  = 
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as  pseudo-cubic  splines  because  of  the  analogy  with  the  cubic 
splines  in  one  variable.  They  belong  to  C^(R),  Pseudo  quintic 
splines  etc.  are  also  considered  in  Duchon  [72,73]. 

A  similar  program  has  been  carried  out  by  Mansfield  [133- 
137]  for  some  spaces  of  smooth  functions  defined  on  a  rectangle 
H.  In  [136]  she  considers  a  space  of  functions 
where  m  and  n  are  positive  integers  and  a  <  CC  <  b,  c<p<d. 
This  space  is  actually  defined  by  completion  of  a  set  of  tensor 
product  functions  with  respect  to  an  appropriate  inner-product, 
and  we  do  not  want  to  define  it  precisely  here.  A  function 
f  e  t'"’ (0!, P)  has  the  following  properties,  however: 


0.25)i 


(i,j) 

e  C(H),  i  = 

0,1, 

« •  •  ^  tn*  1 

and  j  =0,1, 

•  •  •  ^  D' 

(s-j- 

(x,p)  €AC[ 

a,b] 

and  f^^ 

(x,p)  eL^ 

[a,b] 

j  =  0,1,.. 

.,n-l 

(i,s- 

(a,y)  €AC[ 

c,d] 

and  f^^' 

s-i)  ,  V  t2 
(a, y)  e  L 

[c,d] 

i  =  0,1,.. 

.  ,m-l 

(m-1, 

e  AC(H) 

and 

f(m,n)  ^ 

L^(H), 

AC 

stands  for  the 

space  of  absolutely  continuous 

tions  and  where  s  =  m  +  n.  By  constructing  an  appropriate  re¬ 
producing  kernel,  she  is  able  to  solve  problem  (3.19)  with 


(3.26)  e(f)  =  //[f^™'"^  ]^  +  L  f  [f^®'j’ (x,P)  ]^dx 
H  j=0  a 


/e'  /[f^^^^-^>(a,y)]2dy. 

i=0  c 


In  [133],  Mansfield  carries  out  a  similar  analysis  for  a 
space  of  functions  defined  on  the  rectangle  H.  Here 

b]  X  L2[c,d],  where  L2[a,b]  is  the  usual  Sobolev 
space  of  functions  with  absolutely  continuous  derivatives  up  to 
order  m-1,  and  with  f^"*^  e  L^[a, b].  By  constructing  an 


appropriate  reproducing  kernel,  she  now  solves  problem  (3.19) 
with 

(3.27)  e(f)  =  ^  Z  /  [f^"'^j\x,c)]^dx 

H  j  0  a 

m-1  d  .  .  „ 

+  Z  /  [f^^'"^a,y)]"dy  . 

i=0  c 

The  solution  turns  out  to  be  l  piecewise  polynomial  of  degree 
2m- 1  in  x  and  of  degree  2n-..  in  y.  It  is  also  in 
^2m  2,2n  2^^^^  particular  case  of  gridded  data,  it  re¬ 

duces  to  the  tensor  product  of  one-variable  splines  (cf.  the 
following  section) .  Other  more  general  definitions  of  0  are 
also  considered  (with  minor  modifications  on  the  one-dimensional 
integrals) . 

A  more  algebraic  approach  to  constructing  multidimensional 
spline  functions  (which  also  involves  certain  kernel  functions) 
has  been  taken  by  Schaback  [173-174].  His  two-dimensional  ker¬ 
nel  function  is  obtained  as  a  tensor  product  of  one- dimensional 
kernels. 

3.5.  Spline  interpolation  (gridded  data) .  The  problem  of  con¬ 
structing  interpolating  splines  in  two  dimensions  with  gridded 
data  as  in  (3.4) -(3.6)  is,  of  course,  a  special  case  of  the 
general  problems  discussed  in  subsection  3.4.  The  development 
of  the  gridded  data  case  predated  the  more  general  development 
and,  moreover,  is  considerably  simpler.  There  are  a  great  many 
papers  on  two-dimensional  polynomial  splines  and  generaliza¬ 
tions.  We  do  not  have  space  here  to  discuss  all  of  them  in  de¬ 
tail.  We  shall  be  content  to  quote  some  of  the  papers  and  to 
give  a  somewhat  more  complete  discussion  of  polynomial  splines, 
which  are  the  most  widely  used  splines  for  this  problem. 

Some  early  papers  dealing  with  two-dimensional  interpola¬ 
ting  splines  Include  Birkhoff  and  de  Boor  [26],  Birkhoff  and 
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Garabedian  [27],  Price  and  Simonson  [159],  and  Theilheimer  and 
Starkweather  [191],  In  [26]  certain  t icubic  splines  were  intro¬ 
duced  which  were  later  studied  in  detail  in  de  Boor  [32],  The 
problem  was  to  minimize 

(3.28)  /  /  [f^^’^\x,y)]^dxdy 

a  c 


over  all  appropriately  smooth  functions  on  the  rectangle  H 

which  interpolate  the  gridded  data  (3.4)- (3.6).  It  was  found 

that  the  solution  of  this  problem  was  a  certain  bicubic  func- 

2 

tion  with  global  smoothness  C  (H) .  This  problem  was  genera¬ 
lized  to  minimizing 
b  d 

(3.29)  e(f)  =  /  /  [f  (x,y)  ]  dxdy,  m  =  2p,  n  =  2q 
a  c 


in  Ahlberg,  Nilson  and  Walsh  [1,2],  whose  solution  involves 

certain  higher-order  polynomial  splines.  Since  they  are  widely 

used,  we  give  a  short  algebraic  treatment  here. 

The  points  define  a  partition  of  the 

intervals  [a,b]  and  [c,d],  respectively  (cf.  (3.5;).  Suppose 

now  that  x,  <..,<x,<a<b<x,  -  <...<x,  ,  and 

1-m  —  —  -1  —  —  k+2  —  -  K+m-l 


y^.ji  <  •••  5y.ilc<d<  y^^2  5  •  •  '  5  y^+n-1®^®  chosen  arbi¬ 
trarily,  Let  [N?]^  be  the  B-splines  associated  with  the 
1  i-m 

x-partition,  and  let  the  B-splines  associated  with  the  y-parti- 
tion  be  denoted  by  (Nj(y)  For  a  complete  discussion  of 

B-splines  and  their  properties,  see  de  Boor  [36]  in  this  volume 
(or  [33]),  Let 


(3.30)  N^.(x,y)  =  N“(x)N^(y), 

The  linear  space 


(3.31)  ^  =  sp.n 


i  =  1-m,  ...,k  and  j  =  1-n, 


Z 


is  clearly  of  dimension  (k+m) (i+n) .  We  may  now  construct  an 
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element  in  (3.31)  which  interpolates  to  uhe  gridded  data. 

Since  there  are  only  (k+2)(^+2)  data  points  on  the  grid  (cf. 
(3. 4) -(3. 6),  it  is  clear  that  if  we  use  S  to  interpolate,  we 
have 


(3.32)  (k-Hn)(i+n)-(k+2)(i+2)  .  (k+2)  (n-2) +(i+2)  (m-2) +(n-2)  (m-2) 


free  parameters.  Thus,  to  uniquely  define  a  spline,  one  must 
add  additional  conditions.  Recall  that  m = 2p  and  n  =  2q.  Then 
we  might  add  the  extra  conditions 

,y  )  =  y  ).0,  j=0,l,...,m 

(3.33)  ^  J  3  V.p,  ...,m-2 

s^°’^\x^,yQ)  =  (x^,y^^P  =0,  i  =  0, 1, . . . ,k+l 

1-1  =  q,  . . . ,  n~  2 


and 


(3.34) 


I  Cx  V  )  -  3  (x  y  )  -  s  fx  V  ) 

-  vXQ,yQj-s 


=  s 


(V,p) 


^\+l’^£+l^  V-p,...,m-l 

M  -  q, . . . ,  n“  1 . 


These  are  called  the  natural  boundary  conditions,  ami  it  can  be 
shown  that  the  system  of  equations 


(3.35) 


k  i 


Z  Z  a  N  _ 
i=  l-m  j=l-n  ^  ^3 


<vV 


a  =  o,l,  ...,k+l 
P  =  0, 1, , .  .fZ+1 


coupled  with  the  conditions  (3.33) - (3.34)  provides  a  nonsingular 
system  of  equations  for  the  coefficients  This  system 

has  convenient  bandedness  properties  if  the  equations  are  ar¬ 
ranged  properly.  The  resulting  spline  is  precisely  the  solutdbn 
of  the  minimization  problem  (3.29).  The  boundary  conditions 

(3 .33) -  (3.34)  are  the  natural  ones  associated  with  the  problem 
(3.29).  However,  it  is  also  possible  to  specify  lower-order 
derivative  information  along  the  boundary  and  also  obtain  a 
nonsingular  system  of  equations.  The  resulting  spline,  called 
Type  I,  can  also  be  shown  to  satisfy  an  appropriate  minimization 


problem.  However,  for  data- fitting  purposes,  to  use  the  inter- 
polant  with  boundary  derivative  data  one  would  first  have  to 
perform  a  first-stage  approximation  to  find  estimates  for  the 
required  derivatives. 

The  best-known  case  of  the  above  spline  interpolation  is 

the  case  m=n=4,  i.e.,  bicubic  spline  interpolation.  In 

this  case  the  surface  constructed  is  a  piecewise  bicubic  with 

2 

global  smoothness  C  (H)  .  The  natural  boundary  conditions  set 
second-derivative  values  to  0.  Programs  for  computing  natural 
bicubic  interpolating  splines  can  be  found  in  the  IMSL  Library 
[117]  in  FORTRAN.  FORTRAN  programs  for  Type  I  bicubic  splines 
can  be  found  in  Koelling  and  Whitten  [121],  where  the  required 
boundary  information  is  assumed  to  be  input.  An  ALGOL  program 
for  computing  Type  I  bicubic  splines  in  which  boundary  data  are 
automatically  computed  by  fitting  one- dimensional  splines  ap¬ 
pears  in  Spath  [183]. 

Bicubic  spline  interpolation  has  been  widely  applied.  For 
some  references  in  the  Geology  literature,  see  eg.  Anderson  [7], 
Bhattacharyya  [22],  Holroyd  and  Bhattacharyya  [115],  Koelling 
and  Whitten  [121],  and  Whitten  and  Koelling  [206]. 

Problem  (3.29)  has  been  widely  generalized  in  the  spline 
literature.  Instead  of  minimizing  ordinary  derivatives,  one 
may  introduce  general  linear  operators,  and  instead  of  dealing 
with  point  evaluation  functionals,  more  general  linear  function¬ 
als  may  be  permitted.  To  list  some  (but  by  no  means  all)  papers 
dealing  with  such  generalizations,  we  mention  Arthur  [8,9], 
Birkhoff,  Schultz  and  Varga  [29],  de  Boor  [34],  Delvos  [65,66], 
Delvos  and  Schempp  [68,69],  Delvos  and  Schlosser  [70],  Fisher 
and  Jerome  [78,79],  Haussmann  [100],  Haussmann  and  Munch  [104], 
Munteanu  [143,144],  Nielson  [148,150],  Ritter  [164,165],  Sard 
[?71,172],  Schoenberg  [176],  Schultz  [177,178],  Spath  [184,185], 
and  Zavialov  [209-212],  On  L-shaped  regions  and  other  polygons 
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see  Birkhoff  [25]  and  Carlson  and  Hall  [44-49]. 

We  close  this  section  by  mentioning  another  direction  of 
generalization  which  has  led  to  a  considerable  development,  the 
idea  of  spline  blending.  These  methods  are  useful  for  construc¬ 
tion  of  a  surface  which  interpolates  not  only  function  values 
at  isolated  points  but  on  the  grid  lines  themselves;  i.e., 


f  (x,yj) 

=  F(x,yj) 

a  < 

VI 

X 

and  j  =  0,1,.,, , i+1 

f(x.,y) 

=  F(x.,y) 

c  < 

y  <  d 

and  i  =0,1,.. .,k+l. 

To  use  such  blending  methods  one  must  have  F  defined  on  the 
grid  lines.  Thus,  the  methods  could  be  of  value  as  second-stage 
processes.  We  do  not  have  space  to  go  into  detail  on  spline- 
blended  methods.  We  refer  to  the  recent  book  of  Barnhill  and 
Riesenfeld  [20]  for  a  collection  of  papers  on  the  subject  and 
for  further  references.  See  also  the  papers  of  Gordon  [84-87] 
and  Gordon  and  Hall  [88].  Recently,  considerable  effort  has 
gone  into  showing  that  blending  methods  also  arise  as  solutions 
of  appropriate  variational  problems;  see  the  papers  of  Delvos 
[65],  Delvos  and  Rosters  [66],  and  Delvos  and  Malinka  [67]. 

4 .  Local  interpolation  methods 

The  interpolation  methods  discussed  in  section  3  were  glo¬ 
bal  in  nature — that  is,  the  value  f(x,  y)  of  the  constructed 
surface  at  any  given  point  (x, y)  in  D  depends  on  the  values 
of  all  of  the  data  points.  This  generally  means  that  to  compute 
a  representation  for  f  one  has  to  solve  a  fairly  large  system 
of  equations,  and  to  evaluate  f(x, y)  one  generally  has  to 
carry  out  a  considerable  amount  of  arithmetic.  In  this  section 
we  shall  consider  local  schemes  where  the  surface  depends  only 
on  nearby  data  points.  Then  the  construction  will  usually  lead 
to  (a  possibly  large  number)  of  small  systems  of  equations,  and 


moreover,  the  evaluation  of  the  surface  at  a  given  point  will 
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usually  involve  very  little  computation. 

Many  of  the  schemes  mentioned  in  section  3  can  be  made  lo¬ 
cal  in  nature  by  the  following  simple  approach.  Suppose  that 

d 

the  domain  D  is  partitioned  into  subdomains:  D  =  U  D^.  We 
then  seek  a  surface  in  the  form  ^  ^ 

(4.1)  f(x,y)  =  {f^(x,y),  (x,y)  e  D.,  i  =  1,2, ..,,d. 

To  construct  each  individual  f^,  we  suppose  that  are  do¬ 

mains  containing  D^,  which  contain  only  points  which  are  "near" 
D^.  Then  we  use  the  data  (and  only  the  data)  in  to  con¬ 

struct  fj^.  Usually,  we  can  choose  =  D^.  In  most  cases 
the  most  convenient  choices  for  the  subdomains  are  trian¬ 

gles  and  rectangles.  We  discuss  these  two  cases  first. 

4.1.  Triangular  subregions  (scattered  data) .  Suppose  that  we 
are  given  data  at  points  i  =  1,2,  .,.,N  scatter¬ 

ed  throughout  the  plane,  and  let  D  be  the  convex  hull  of 
these  points.  It  is  more  or  less  clear  that  by  drawing  lines 
from  point  to  point  we  can  construct  a  set  of  triangles  with 
vertices  at  the  which  partition  D.  It  is  also  clear  that 

given  any  set  of  points,  this  triangularlzatlon  of  D  is  not 
usually  uniquely  defined  (see  Figure  2  below  for  two  different 
triangularizations  of  the  same  region) .  Moreover,  as  the  fig¬ 
ure  shows,  some  triangularizations  are  superior  to  others  in 
the  sense  that  they  exhibit  fewer  of  the  less  desirable  long 
thin  triangles. 

I 

i 

I  _  Figure  2.  Triangularlzatlon 
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The  design  of  an  algorithm  to  divide  a  region  into  accept¬ 
able  triangles  with  vertices  at  prescribed  points  is  not  as 
easy  as  it  sounds.  Two  algorithms  in  the  literature  which  are 
designed  to  give  good  triangularizations  can  be  found  in  Caven¬ 
dish  [50]  and  in  Lawson  [128]. 

The  simplest  approach  to  defining  a  local  interpolating 
surface  is  to  construct  f^(x,y)  to  be  of  the  form  a^^  +  a2X  + 
a^y  in  each  triangle.  The  data  at  the  three  comers  of  the 
triangle  determine  the  coefficients  for  that  piece  of  f  (the 
corresponding  system  will  be  nonsingular  provided  the  triangle 
is  nondegenerate) ,  This  procedure  produces  a  piecewise  linear 
surface  which,  in  fact,  will  be  globally  continuous.  This  last 
property  follows  from  the  fact  that  along  the  sides  of  the  tri¬ 
angle  the  functions  reduce  to  straight  lines  joining  the  ver¬ 
tices.  This  method  has  been  used  by  several  authors  for  data 
fitting,  e.g.,  Lawson  [128]  and  Whitten  [206].  For  some  con¬ 
touring  routines  based  on  this  local  interpolation  scheme,  see 
section  8. 

If  we  desire  to  interpolate  several  sets  of  data  defined 
on  the  same  triangularization,  it  may  be  more  convenient  to 
compute  Lagrangian  functions  rather  than  to  compute  the  surface 
in  each  triangle  separately.  In  particular,  it  is  clear  that 
we  can  construct  functions  {0j(x,  y))^  with  the  property 

(4.2)  0j(Xj^,y^)  =  i,j  =  1,2,  ...,N. 

These  functions  can  be  constructed  as  pyramids  in  such  a 
way  that  the  function  0.  has  support  only  on  the  triangles 
surrounding  the  point  (Xj,yj)  (see  Figure  3).  In  terms  of 
these  Lagrangian  functions,  the  interpolating  surface  is  given 
by 

N 

(4.3)  f(x,y)  =  Ef0  (x,y). 

j=l  J  J 


Figure  3.  A  Lagrange  Element 

The  Lagranglan  approach  to  local  interpolation  is  very 
reminiscent  of  the  finite  element  method  in  which  the  solution 
of  an  operator  equation  is  sought  in  the  form  of  a  linear  com¬ 
bination  of  a  set  of  functions  (called  elements)  with  the  pro¬ 
perty  (4.2).  (See  e.g.,  Prenter  {157],  Schultz  [179],  or 
Strang  and  Fix  [188].)  There  is  no  need  to  restrict  the  ele¬ 
ments  to  be  piecewise  linear  functions — we  may  use  higher-order 
polynomials,  rational  functions,  or  even  more  complicated  func¬ 
tions.  In  fact,  if  we  are  careful  in  the  construction,  we  may 
be  able  to  construct  elements  with  small  support  but  higher 
global  smoothness. 

There  are  a  great  many  papers  in  the  finite-element  litera¬ 
ture  concerned  with  defining  convenient  smooth  elements  (La- 
grangian  functions  with  small  support) ,  To  mention  a  few,  see 
Barnhill,  Birkhoff,  and  Gordon  [16],  Barnhill  and  Gregory  [17, 

18],  Barnhill  and  Mansfield  [19],  Birkhoff  and  Mansfield  [28], 
Bramble  and  Zlamal  [39],  Goel  [83],  Hall  [94],  Mitchell  [141], 
Mitchell  and  Phillips  [142],  Nicolaidis  [146,147],  Zenisek  [213], 
Zienkowicz  [214],  and  Zlamal  [215-217].  The  books  on  finite 
elements  of  Aziz  [13],  de  Boor  [35],  Strang  and  Fix  [188],  and 
Whiteman  [198]  should  also  be  consulted. 
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The  construction  of  elements  with  higher-order  smoothness 
becomes  increasingly  difficult.  For  example,  it  is  shown  in 
Mansfield  [137]  that  to  get  an  element  with  support  on  the  tri¬ 
angles  surrounding  and  with  global  continuity  C  (D) ,  it 

is  necessary  to  use  polynomials  of  degree  5  at  ’east.  (Matters 
are  somewhat  simpler  on  regular  triangularizations,  see  subsec¬ 
tion  4.2  below.) 

We  close  this  subsection  by  mentioning  that  it  is  also  pos¬ 
sible  to  perform  interpolation  using  elements  based  on  triangles 
to  data  which  also  involves  derivatives,  or  in  analogy  with  the 
blending  methods,  to  data  which  includes  function  values  along 
the  edges  of  the  triangles.  (See  e.g.,  Barnhill,  Birkhoff,  and 
Gordon  [16],  or  Barnhill  and  Gregory  [17,18],)  These  methods 
are  not  directly  applicable  to  the  scattered  data  Problem  1.1, 
but  may  be  useful  as  second-stage  methods. 

4.2.  Regular  triangularizations.  When  the  data  is  distributed 
such  that  the  region  can  be  triangulated  into  a  set  of  congru¬ 
ent  triangles,  then  it  Is  extremely  advantageous  to  use  the  La¬ 
grange  approach.  In  particular,  in  this  case  we  can  find  an 
element  0  with  value  1  at  (0,0)  such  that  all  other  elements 
are  translates  of  0.  In  this  case,  f  takes  the  form 
N 

(4.4)  f(x,y)  =  E  F  0((x,y)  -  (x  ,y  ))  . 

j=l  J  J  J 

We  illustrate  this  with  a  couple  of  examples.  Suppose  that 
we  are  given  data  at  points  chosen  from  the  collection 

(4.5)  =  t(i,  j)  ^  ^i,  jeZ  "  ^  "  [integers). 

These  points  lie  on  the  comers  of  a  triangular  grid  as  shown 
in  Figure  4. 

It  Is  shown  in  Zwart  [218, p.  673]  that  there  exists  a  func- 
1  2 

tion  0  e  C  (R  )  which  is  1  at  the  origin  and  0  at  all  other 
points  in  Q,  and  has  support  on  the  shaded  region  in  Figure  4.  ■ 


Figure  4.  A  Regular  Triangularization 


This  function  is  constructed  as  a  piecewise  quadratic  polyno¬ 
mial.  A  similar  element  has  been  constructed  by  Powell  [156] 
(the  figure  on  page  267  of  [156]  should  be  rotated  45°  to  see 
this) . 

To  give  another  example,  suppose  that  we  consider  the  set 
of  points  which  lie  at  the  vertices  of  the  grid  defined  by 

equilateral  triangles  shown  in  Figure  5. 


Figure  5.  Another  Regular  Triangularization 


It  is  shown  in  Fredrickson  [81]  that  there  exists  a  function  0 

which  has  value  1  at  the  origin  and  value  0  at  all  other  points 

2  2 

in  The  function  0  is  in  C  (R  ),  consists  of  piecewise 

quartics,  and  has  support  in  the  region  shown  in  Figure  5, 

Fredrickson  also  constructs  a  piecewise  cubic  element  with  the 

1  2 

same  support  but  which  is  only  C  (R  )  .  For  right  triangles 
see  Carlson  and  Hall  [44]. 


4.3.  Rectangular  subregions.  In  this  section  we  suppose  that 
we  have  data  given  at  points  lying  on  a  rectangular  grid  as  in 
(3. 4) -(3. 6),  and  consider  local  interpolation  methods.  The 
simplest  approach  here  (cf.  the  triangularization  case)  is  to 
construct  a  separate  bilinear  function  f(x, y)  =  a^^  +  a^x  + 
a^y  +  a^xy  in  each  subrectangle,  x 

using  the  four  comer  values  to  determine  the  coefficients. 

Since  the  bilinear  patches  reduce  to  linear  functions  on  the 
grid  lines,  the  global  surface  is  C(R) . 

Several  authors  have  considered  constructing  functions  on 
each  of  the  using  higher-order  polynomials.  This  requires 

additional  information  in  addition  to  the  four  corner  values. 

For  example,  if  one  seeks  a  bicubic 


3  3 

(4.6)  f(x,y)  =  E  Z  a  x\^, 
i=0  j=0 


there  are  16  coefficients  to  determine.  These  could  be  deter¬ 
mined  by  the  four  comer  values,  plus  the  values  of  f  ,  f  , 

X  y 

and  f^^  at  each  comer.  To  determine  these,  one  must  perform 
some  first- stage  process.  For  some  approaches  to  this,  see 
Akima  [5],  Hessing,  et  al  [114],  and  Shu,  et  al  [181].  A  FOR¬ 
TRAN  program  for  Akima' s  method  can  be  found  in  [6],  Nonpoly¬ 
nomial  patches  have  also  been  considered;  e.g.,  see  Birkhoff 
and  Garabedian  [27]. 

The  Lagrange  (finite  element)  approach  can  also  be  used  in 
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the  case  of  rectangular  gridded  data.  In  particular,  if  we  can 
construct  a  function  satisfying  (4,2)  with  local  support,  then 
the  surface  f  given  by  (4.3)  will  interpolate  and  the  method 
will  be  local  in  character.  As  before,  the  Lagrange  approach 
is  especially  convenient  if  the  grid  is  regular,  i.e.,  if  all 
subrectangles  are  congruent.  To  illustrate  this,  suppose 

that  the  are  actually  the  unit  squares;  i.e.,  the  data 

points  lie  in  the  set 


(4.7)  =  ((i,j))  i,j  e  Z,  Z  =  {integers}. 

T'r  get  a  quadratic  element,  we  may  simply  rotate  the  ele¬ 

ment  of  Zwart  [218]  considered  in  the  last  section  by  45  degrees 
(cf.  Figure  4),  or  we  may  take  the  element  of  Powell  [156]. 


4.4.  Parametric  representations.  The  methods  discussed  in  the 
last  section  is  concerned  with  data  given  on  a  rectangular  grid. 
By  using  parametric  representations,  it  is  possible  to  construct 
similar  local  interpolating  surfaces  for  data  given  at  the  cor¬ 
ners  of  any  partition  of  D  consisting  of  quadrilaterals.  In 
this  section  we  briefly  describe  how  this  might  proceed. 

Suppose  Q  is  a  particular  quadrilateral  subregion  of  D 
of  interest.  In  addition,  suppose  that  x(s,t),  y(s,  t),  and 
z(s,t)  are  functions  defined  on  the  unit  square  U  =  [0,1]  x  [0,1] 
with  the  properties  that  as  (s,  t)  runs  over  the  boundary  of 
U,  (x(s,  t)  ,  y  (s,  t)  )  runs  over  the  boundary  of  the  quadrilateral; 
the  four  comers  of  U  correspond  to  the  four  comers  of  Q; 
and  z(s,  t)  takes  on  the  desire  J  data  values  at  the  four  cor¬ 
ners  of  U.  In  this  case,  the  triple  (x(s,  t) ,  y  (s,  t)  ,z  (s,  t)) 
provides  a  parametric  representation  of  a  piece  of  surface  de¬ 
fined  over  Q  interpolating  the  data. 

The  problem  of  constructing  parametric  representations  of 
interpolating  functions  has  been  considered  in  a  number  of  pa¬ 
pers.  Several  papers  on  these  methods  and  a  host  of  references 
can  be  found  in  the  book  of  Barnhill  and  Riesenfeld  [20];  see 


also  the  survey  paper  of  Shu  et  al  [181],  Such  surfaces  are 
sometimes  called  Coon's  surfaces,  cf.  Coons  [59],  and  are  of 
considerable  interest  in  the  field  of  computer-aided  geometric 
design.  To  mention  just  a  few  of  the  actual  papers,  see  Ahuja 
and  Coons  [4],  Eamshaw  and  Youille  [74],  Ferguson  [77],  Hayes 
[107],  Hosaka  [116],  and  Mangeron  [132], 

There  also  has  been  some  effort  directed  towards  construct 
ing  elements  (Lagrange  functions)  associated  with  other  less 
regular  subsets  of  the  plane.  We  mention,  for  example,  the 
work  of  Ciarlet  and  Raviart  [55],  Wachspress  [194,195],  and 
Zlamal  [217]  in  which  elements  are  constructed  for  domains 
involving  curved  edges. 


4.5,  Local  Shepard  methods.  It  is  possible  to  modify  the  meth 
od  discussed  In  subsection  3.3  to  make  it  local.  For  example, 
following  Shepard  [180],  suppose  we  fix  0  <  R  and  define 


r 


l/r 


0  <  r  <  f , 


(4.8)  xKr)  =  { 


II  (|  -  1)  R/3  <  r  <  R, 

0  ,  R  <  r  . 


This  function  is  continuously  differentiable  and  vanishes  iden- 

R, 


tically  for  r  <  R,  Now  with  r^  as  in  (3.8),  we  define 


(4.9)  f(x,y)  =  < 


N 

2  F  [\l/(r.) 

1=1  ^  ^ 

N 

Z  [^^(r.)]^ 

i=l  ^ 


,  when  r^  ^  0,  all  i 


,  when  r  =  0 . 


2 

Formula  (4.9)  is  defined  at  all  (x, y)  in  the  plane  R  . 
By  definition  it  interpolates  the  values  F^  at  the  data 
points  (Xj^,y^),  i  =  1,2,  ...,N,  The  values  at  non-data  points 
are  obtained  as  weighted  averages  of  the  data  values  F^,  but 
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I; 


■A'' 


only  those  which  lie  at  points  within  a  distance  of  R  of 
(x,y).  Thus,  the  formula  is  local. 

To  use  this  method  in  practice  it  is  necessary  to  choose 
a  reasonable  value  for  R.  The  aim  is  to  find  R  so  that  for 
every  (x,y)  a  reasonable  number  ol  data  points  will  fall  in 
the  disk  centered  at  (x, y)  of  radius  R  .  It  would  also  be 
possible  to  let  R  depend  on  (x, y) ,  i.e.,  to  use  different 
valuffi  of  R  in  different  subregions  of  D. 

5 .  Global  approximation 

As  mentioned  in  the  introduction,  frequently  the  data  does 
not  warrant  constructing  an  interpolating  function  (e.g.,  be¬ 
cause  of  errors) .  In  such  cases  it  may  be  preferable  to  con¬ 
struct  a  surface  which  only  approximates  the  data.  In  this  sec¬ 
tion  we  discuss  some  global  approximation  methods. 


5.1.  Polynomial  least  squares.  The  general  theory  of  discrete 
least-squares  fitting  is  very  well  known.  To  briefly  review, 
suppose  that  ^  given  functions  on  D.  Define 


(5.1)  *(a)  = 


N  n 

Z  Z  a  .0  .  (x  .y.)  -  F. 
.  1  .  ,  1  1  i  1  1 

1-1  j=i 


T  n 

(a^, ...,a^)  is  any  vector  in  R  . 


where  a 
lem  is  to  find  a*  such  that 


Then  the  prob- 


(5.2)  $  (a*)  =  min  <I>  (a)  . 

a 

The  corresponding  function 
n 

(5.3)  f(x,y)  =  Z  a*  0  (x,y) 

j=l  J  J 

is  called  the  discrete  least-squares  approximation  of  the  data 
{Ff)^.  Usually  one  takes  n  considerably  smaller  than  N.  In 
this  section  we  briefly  discuss  least  squares  using  polynomials. 
Before  doing  so,  however,  we  make  a  few  general  remarks  about 
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are  orthonormal  with 


solving  the  general  least-squares  problem. 

There  are  several  approaches  to  solving  (5,2).  Perhaps 
the  neatest  is  the  case  where  the 
respect  to  the  inner-product 
N 

(5.4)  ((^,\lr)  =  Z  (»(x  ,y  )t  (x  ,y  ) . 

i=l  ^ 


Then  the  solution  of  (5.2)  can  be  written  down  explicitly  as 
n 

(5.5)  f (x,y)  =  Z  F  0  (x,y)  . 

j=l  ^  ^ 

A  second  very  well-known  approach  to  solving  (5.2)  is  via 
the  normal  equations 

(5 . 6)  A*A  a  =  A*F  , 

T 

where  F  =  (Fj^,  ...,F^)  is  the  vector  of  data  values,  and  where 

(5.7)  A  =  (0j(x^,yp) 

In  some  cases  the  normal  equations  are  a  perfectly  acceptable 
way  to  compute  least- squares  approximation,  but  in  other  cases 
the  system  (5.6)  may  be  ill-conditioned  (or  even  singular--cf . 
the  following  subsection  for  spline  least  squares) .  This  ap¬ 
proach  is  also  not  convenient  should  side  conditions  be  desired 
(e.g.,  by  imposing  actual  interpolation  at  some  of  the  values). 
For  more  on  the  normal  equations,  see  any  book  on  Numerical 
Analysis, 

A  more  modem  method  of  solving  least-squares  problems  is 
to  use  general  matrix  methods.  Specifically,  consider  the  ob¬ 
servation  equations 

(5.8)  Aa  =  F. 

It  can  be  shown  that  by  applying  a  series  of  matrix  transfor¬ 
mations  to  this  system,  one  can  obtain  a  set  of  equations  giving 
the  vector  a*.  For  a  complete  description  of  methods  of  this 


type  see  Lawson  and  Hanson  [129j  or  Stewart  [187],  Matrix 
methods  are  quite  amenable  to  the  adding  of  side  conditions  and 
can  also  be  designed  to  take  account  of  rank-deficiency  of  the 
matrix  A  (which  corresponds  to  the  case  of  singular  normal 
equations) . 

Polynomial  discrete  least-squares  fitting  has  been  widely 
used  for  fitting  surfaces  to  data,  both  scattered  and  regular. 
Several  authors  have  developed  algorithms  for  polynomial  dis¬ 
crete  least-squares  fitting  cf  scattered  data  by  constructing 
orthonormal  polynomials  (e.g,  by  Gram-Schmidt  orthonormaliza¬ 
tion)  .  See,  for  example,  Cadwell  and  Williams  [42],  Crain  and 
Bhattacharyya  [61],  and  Whitten  [201,202],  The  latter  contains 
a  FORTRAN  program. 

When  the  data  are  more  regularly  distributed,  polynomial 
least-squares  fitting  can  often  be  simplified.  For  example,  if 
the  data  lie  on  a  grid  as  in  (3. 4) -(3, 6),  then  the  desired  or¬ 
thogonal  polynomials  are  simply  products  of  the  one-dimensional 
orthogonal  polynomials  associated  with  the  one-dimensional  inner 
products  corresponding  to  3nd  respectively,'  e,g,, 

see  Cadwell  [41]  or  Clenshaw  and  Hayes  [56],  as  well  as  the  sur¬ 
vey  papers  of  Hayes  [105,108,109], 

There  are  also  special  methods  for  handling  data  which  are 
not  on  a  grid  but  instead  lie  on  parallel  straight  lines.  For 
example,  Clenshaw  and  Hayes  [56]  have  developed  methods  using 
expansions  in  terms  of  Tchebycheff  polynomials  (although  the 
method  actually  only  produces  an  approximation  to  the  least- 
squares  fit  rather  than  the  actual  minimum) , 

Polynomial  least  squares  can  also  be  interpreted  as  multi¬ 
dimensional  regression  as  practiced  by  statisticians,  e,  ,, 
see  Effroymson  [75],  For  example,  if  we  are  trying  to  fit  a 
function  in  the  form 
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dx  dy  .  . 

f (x,y)  E  Z  a  X  y-*, 
i-0  j-0  -■ 

then  by  defining  new  variables  by 

V  H  V  =  0, 1, . ,  .,dx 

""vCdy+D+M  -  y  >  .  0,1,... ,dy 

we  can  write 

d 

y)  =  E  b  z  d  =  dxdy  +  dx  +  Iv, 

1=0  ^  ^ 

and  the  problem  becomes  one  of  fitting  a  linear  function  in 
several  variables. 

We  close  this  section  by  observing  that  in  some  cases  it 
may  be  desirable  to  consider  weighted  least  squares.  In  parti¬ 
cular,  if  we  have  positive  weights  w^  >0,  i  =  1, 2,  ...,N,  then 
we  may  replace  4'  in  (5 . 1)  by 

It  is  interesting  to  note  that  the  interpolation  formula 
of  Shepard  discussed  in  section  3.3  can  be  interpreted  in  terms 
of  weighted  least-squares  fitting.  In  particular,  fix  (x, y) 
in  D,  and  let  r^(x, y)  be  the  distance  from  (x, y)  to  the 
point  (x^,yj^)  as  before.  Now  set  w^  =  r^*^,  and  consider 
least- squares  approximation  by  a  constant  c,  using  these 
weights.  Then  one  easily  computes  that  the  least-squares  choice 
of  c  is 


N 

N 

E  w,F, 

E  F,  rT' 

1  ‘  ‘ 

1  ^  " 

N 

N 

E  w. 

E  r^ 

1  i 

1  i 

This  approach  was  adopted  by  Pelto,  Elkina  and  Boyd  [152]  (as 
pointed  out  to  me  by  Chuck  Duris) . 


5.2.  Discrete  least-squares  fitting  by  splines.  As  outlined 
in  the  previous  subsection,  discrete  least  squares  can  be  car¬ 
ried  out  with  any  finite  set  of  functions.  It  is  not  surpris¬ 
ing  that  a  number  of  authors  have  tried  using  tensor  product 
splines.  See,  e.g.,  Halliday,  Wall,  and  Joyner  [96],  Hayes  and 
Halliday  [110],  Jordan  [119],  Hanson,  Radbill,  and  Lawson  [97], 
and  Whiten  [199].  Hayes  and  Halliday  have  developed  both  ALGOL 
and  FORTRAN  programs.  It  is,  on  the  other  hand,  perhaps  some¬ 
what  surprising  that  least-squares  fitting  with  splines  can  be 
somewhat  problematical.  We  briefly  discuss  the  method. 

Suppose  that  H  =  [a, b]  x  [c,d]  is  a  rectangle  containing 
the  domain  D  of  interest.  Let 


and 


{yj}Q^^  be  parti- 
,k+l  i+1 


tions  of  [a,b]  and  [c,d],  respectively,  and  let  , 

1 j  l-m, i-n 

be  the  tensor  product  B-splines  discussed  in  section  3.5.  We 
consider  discrete  least- squares  fitting  using  these  (k+m)  (i+n) 
B-splines. 

To  explain  how  difficulties  can  arise  with  spline  least- 
square  fitting,  we  observe  that  it  is  quite  easy  for  the  matrix 
A  in  the  observational  equations  (5.8)  to  be  rank-deficient. 

On  a  trivial  level  this  can  happen  if  for  some  B- spline 
none  of  the  data  points  lies  in  its  support.  This  deficiency 
can,  of  course,  be  easily  removed  by  dropping  this  particular 
B-spline  from  the  set  being  used  to  approximate.  But  rank  de¬ 
ficiency  can  also  occur  in  more  subtle  ways  because  of  the 
local  support  properties  of  the  functions.  This  problem  can  be 
overcome  with  properly  designed  algorithms.  See  Hayes  and  Halli¬ 
day  [110]  for  a  careful  discussion  of  spline  least-squares  fit¬ 
ting.  Lawson  and  Hanson  [129]  include  a  general  discussion  of 
how  to  handle  rank  deficient  least-squares  problems. 

If  we  operate  in  terms  of  the  normal  equations,  then  it 
may  well  occur  that  the  normal  equations  are  in  fact  singular. 
This  is  again  due  to  the  local  property  of  the  B-splines  com- 
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blned  with  the  discrete  inner- product.  Even  when  it  is  not 
singular,  the  set  of  normal  equations  can  be  ill-conditioned 
(even  though  it  is  a  relatively  sparse  matrix  with  a  kind  of 
repeated  band- s true tu re)  . 

Discrete  least  squares  can  also  be  carried  out  with  vari¬ 
ous  finite  dimensional  linear  spaces  of  blended  functions.  For 
an  extensive  study  of  such  methods,  see  the  dissertation  of 
Doty  [71]. 

5.3.  Discrete  and  It  approximation.  Instead  of  performing 
discrete  least  squares,  we  may  consider  the  following  discrete 
approximation  problem:  Given  functions  (0.)^  defined  on  D, 
we  seek  a*  so  that 

N  n 

(5.9)  4*(a)  =  Z  I  E  a  0  (x  ,y.)  -  F  | 

i=l  j=l  J  J  ^  ^  ^ 

is  minimized.  Alternatively,  we  may  minimize 

n 

(5.10)  'I>(a)  =  max  |  E  a  0  (x  ,y  )  -  F  |, 

lSi§N  j=l  J  J  ^  ^  ^ 

These  are  the  usual  and  best  approximation  problems. 

Both  of  these  problems  can  easily  be  reformulated  as  linear 
programming  problems  for  the  determinations  of  the  optimal  a* 
(cf.  Rabinowicz  [160,161]  or  Rosen  [167]).  Reasonable  choices 
for  the  {0j]  would  be  low-degree  polynomials  if  D  is  small, 
or  possibly  spline  functions. 

Discrete  approximation  methods  of  this  type  have  had  rela¬ 
tively  little  exposure  in  the  literature.  For  some  results 

using  tensor  product  splines  in  the  Z  problem,  see  Rosen. 

00 

The  optimal  a*  was  obtained  there  by  using  the  standard  sim¬ 
plex  method  on  the  associated  dual  linear  programming  problem. 

The  Z  problem  can  also  be  solved  by  using  Remez-type 

00 

algorithms.  For  an  algorithm  which  performs  generalized 


rational  approximation  (and  thus  can  also  be  used  for  polynomi¬ 
al  approximations)  see  Kaufman  and  Taylor  [120].  Theoretical 
considerations  for  Tchebycheff  approximation  in  several  vari¬ 
ables  can  be  found  in  Collatz  [58]  or  Weinstein  [196],  for  ex¬ 
ample. 

5.4.  Spline  smoothing  (scattered  data) .  In  this  section  we 
consider  some  minimization  problems  similar  to  those  discussed 
in  section  3.4,  but  where  the  class  of  admissible  functions  is 
not  required  to  interpolate  and  where  the  functional  to  be  mini¬ 
mized  includes  a  term  measuring  how  close  the  function  comes  to 
fitting  the  data.  To  be  more  specific,  suppose  X  is  a  linear 
space  of  "smooth"  functions  and  that  0  is  a  functional  on  X 
which  measures  the  smoothness  of  an  element  in  X.  Suppose  in 
addition  that  E  is  a  functional  defined  on  X  which  measu.:es 
how  well  a  function  fits  the  data.  Then  the  spline- smoothing 
problem  is  the  following; 

(5.11)  Find  s  e  X  such  that  p(s)  =  inf  p(u), 

ueX 

where 

(5.12)  p(f)  =  e(f)  +  E(f). 

The  abstract  theory  of  spline  smoothing  has  been  well 
developed^  see,  e.g.,  the  book  of  Laurent  [127]  and  references 
therein.  To  illustrate  the  ideas,  we  briefly  discuss  a  couple 
of  examples.  We  suppose  as  in  section  3.4  that  X  is  a  semi- 
Hilbert  space  and  that  0  is  a  seminorm  on  X  with  N  = 

{f  e  X;  0(f)  =0}.  We  also  suppose  that  X  is  actually  a 
function  space  defined  on  a  domain  D,  and  that  the  point  eval¬ 
uators  at  {(x^,y^)}^  are  bounded  linear  functionals  on  X. 

We  define 

(5.13)  E(f)  =  p  Z  [f(x  ,y  )  -F  ]% 

i=l  ^  ^  ^ 

where  p  is  a  fixed  positive  constant.  Then  it  can  be  shown 
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(cf.  Duchon  [72,73])  that  the  solution  of  Problem  (5.11)  is  a 
spline  which  can  be  written  in  the  form  (3.20),  where  now  the 
coefficients  are  determined  from  the  linear  system 
N  d 

^  K((Xj,yj)  ;  (x^,y^))a^  +  E  b^P^(Xj,y^)  +aj/p  =  Fj, 

(5.14)  j  =  1,2,...,N, 
d 

E  a  p  (x  ,  y.)  =  0,  k  =  1,2,  ...,d. 

i^l  i  k  1  1 

As  in  section  3.4,  the  application  of  this  method  depends 
on  constructing  a  reproducing  kernel  K.  If  0  is  chosen  as 
in  (3.22),  Atteia  [10-12]  and  Thomann  [192,193]  considered 
spline  smoothing  for  spaces  of  smooth  functions  on  the  rectangle 
and  on  the  disc  (the  latter  even  contains  ALGOL  programs) . 

2 

Duchon  [72,73]  considers  similar  problems  defined  on  D  =  R  . 

A  similar  spline- smoothing  problem  has  also  been  consider¬ 
ed  by  Pivorarova  [154],  where  9  is  taken  to  be 

(5.15)  e(f)  =  //[D^f]^  +  [D^f]^. 

See  also  Kubik  [123]. 

5.5.  Smoothing  splines  (gridded  data) .  In  section  3.5  we  con¬ 
sidered  several  minimization  problems  whose  solutions  led  to 
interpolating  polynomial  splines  (and  generalizations)  .  In  con¬ 
junction  with  the  development  of  interpolating  splines  for 
gridded  data,  there  was  a  concurrent  development  of  smoothing 
splines.  For  example,  instead  of  minimizing  the  integral  0 
in  (3.29)  over  appropriate  smooth  interpolating  functions,  we 
may  minimize  instead  p(f)  =  0(f)  +  pE(f),  where  E  is  given 
by 

k+1  ^+1  n 

(5.16)  E(f)  =  E  E  [f(x  ,y  )  -  F  ]^ 

i=0  j=0  ^  J 

For  results  in  this  direction,  see  e.g.  Nielson  [149,150].  For 
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0  given  by  (3.29),  the  smoothing  splines  are  again  polynomial 
splines.  Again,  more  general  linear  differential  operators  and 
more  general  linear  functionals  can  be  considered. 

5.6.  Continuous  least  squares.  The  method  of  continuous  least 
squares  is  not  directly  suited  to  fitting  surfaces  to  discrete 
data,  but  it  can  be  of  use  as  a  second- stage  process,  so  we 
briefly  review  it.  We  suppose  now  that  F  is  a  function  de¬ 
fined  on  D  which  we  wish  to  approximate,  and  that 
given  functions  on  D,  We  define 

(5.17)  (f,g)  =  //  f (x,y)g(x,y)dxdy,  l|f||^  =  (f,f) 

D 

and 

n  2 

(5.18)  $(a)  =  II  E  a.0.  -  f||  . 

j  =  l  J  J 

The  problem  is  to  find  a*  to  minimize  ^(a).  The  solution  is 
given  by  solving  the  normal  equations 

(5.19)  Aa  =  r, 
where 

A  =  and  r  =  [(0^,F),...,  (0^,F)]^. 

For  reasonably  nice  approximating  functions  it  is  often 
possible  to  compute  the  normal  matrix  exactly.  In  practice, 
the  difficulty  lies  in  evaluating  the  right-hand  sides.  Gener¬ 
ally  a  quadrature  formula  is  required  for  this.  One  advantage 
of  the  method  would  be  that  if  several  data-fitting  problems 
are  to  be  solved  using  the  same  set  of  approximating  functions, 
one  can  do  the  work  of  inverting  the  normal  matrix  just  once 
and  re-use  the  result  as  many  times  as  desired. 

Reasonable  choices  for  the  approximating  functions  Include 
polynomials,  or  better  yet,  tensor  product  B-spllnes  as  in 
(3.30).  Here  the  singularity  problems  do  not  crop  up  for  the 
splines  because  we  are  integrating  instead  of  summing  over 
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finitely  many  points.  The  normal  matrix  in  this  case  has  a  kind 
of  repeated  band  structure.  The  entries  can  be  computed  exactly, 
t)y  Gaussian  quadrature  (cf.  de  Boor,  Lyche  and  Schumaker 
[38]) .  Uniform  best  approximation  by  tensor  products  of  splines 
has  also  been  considered,  e.g.,  see  Sommer  [182]. 

6.  Local  approximation  methods 

As  pointed  out  at  the  beginning  of  se(  :ion  4,  there  are 
many  advantages  which  accrue  if  one  uses  local  methods  rather 
than  global  ones.  In  this  section  we  discuss  some  local  approxi¬ 
mation  schemes. 

6.1.  Patch  methods.  As  in  the  case  of  interpolation,  the  sim¬ 
plest  approach  to  obtaining  local  approximation  methods  is  to 
partition  the  domain  and  to  define  a  surface  (patch)  on  each 
subdomain  separately.  In  particular,  suppose  that  D  =  U{D^)p 
where  are  disjoint  subsets  of  D.  Then  we  may  seek  f  in 

the  form 

(6.1)  f(x,y)  =  {f^(x,y),  (x,y)  e  D^,  i  =  l,2,...,d. 

To  construct  the  patch  f^(x,y),  we  might  use  the  data  available 
in  the  subregion  D^.  In  certain  cases,  however,  it  may  well  oc¬ 
cur  that  no  data  at  all  are  available  in  the  set  D^.  In  this 
case  we  may  choose  a  somewhat  larger  set  of  points  "near" 

D^,  and  use  the  data  in  to  construct  f^^.  For  any  given 

method,  it  should  be  possible  to  make  the  choice  of  adaptive 

so  that  the  size  of  is  kept  as  small  as  possible  consistent 

with  the  amount  of  data  desired  for  the  construction  of  f^. 

The  patch  method  outlined  above  can  be  used  with  any  of  the 
approximation  methods  discussed  in  section  5.  For  example,  one 
might  choose  to  use  polynomials  (of  low  order),  and  to  do  dis¬ 
crete  least-squares  approximation.  Or,  one  might  use  or 

£  approximation  or  some  other  convenient  space  (e.g.  splines) 

00 
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instead  of  polynomials.  The  main  point  is  to  keep  the  size  of 
each  individual  patch  problem  (ind  thus  the  size  of  the  corre¬ 
sponding  system  of  equations)  small.  We  may  have  to  solve  a 
lot  of  systems  of  equations,  but  each  will  be  small  and  fairly 
we 11- conditioned. 

To  illustrate  how  the  adaptive  feature  might  be  implemented, 
suppose  that  the  domain  D  of  interest  has  been  enclosed  in  a 
rectangle  H  and  that  a  partition  of  H  is  defined  by  H  = 


(6.2)  a  =  x„  <  <...<  Xj^^j  =  L,  c  =  yo<yi<--'<yi+i 


d. 


Now  suppose  that  we  want  to  do  discrete  least-squares  fitting 

a  +  bx  +  cy  on  H 


using  a  patch  of  the  form  f^^Cx,  y) 


ij’ 


In  this  case  it  would  be  reasonable  to  require  that  at  least 


3  pieces  of  data  should  be  used  to  construct  f 
does  not  contain  3  pieces  of  data,  we  e;vpand  H 


ij' 

ij 


“  ."ij 

to  by 


adding  all  bordering  rectangles.  If  this  does  not  contain 
3  pieces,  we  again  add  all  bordering  rectangles,  etc.  We  then 
compute  the  discrete  least-squares  polynomial  using  the  data  in 
but  then  we  use  the  resulting  function  only  in  The 

process  may  be  repeated  to  define  each  required  patch.  This 
kind  of  adaptive  algoritlm  is  very  easy  to  program. 

In  using  patch  me  ,i  ods  Lo  get  local  inteirpolatlon  methods, 
we  concentrated  on  methods  using  data  at  comers  of  triangles 
or  rectangles,  and  by  choosing  appropriate  forms  for  the  patches, 
it  was  possible  to  get  the  individual  patches  to  match  together 
to  give  a  continuous  global  surface  (or  with  more  sophisticated 
patches,  even  C^(D)  or  higher).  Here,  however,  where  the  in¬ 
dividual  patches  are  determined  by  approximation,  it  is  not 
very  likely  that  the  patches  will  match  up,  and  the  global  sur¬ 
face  will  generally  not  even  be  continuous.  For  most  applica¬ 
tions,  this  is  a  serious  drawback.  However,  as  we  shall  see  in 


V! 
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section  If  patch  approximation  methods  can  still  be  very  useful 
as  first-stage  methods. 

6.2.  Direct  local  methods.  In  this  section  we  discuss  some 
local  methods  in  which  the  approximating  surface  is  constructed 
directly  from  the  data  without  solving  any  systems  of  equations. 
It  will  be  convenient  to  pose  a  more  general  problem  than  pre¬ 
viously  considered. 

Let  9  be  a  linear  space  of  functions  defined  on  D,  and 
suppose  that  linear  functionals  defined  on  y.  Let 

N 

^^1^1  ^  prescribed  set  of  functions  defined  on  D.  Then  we 

are  interested  in  approximation  schemes  of  the  following  form; 

N 

(6.3)  QF(x,y)  =  ZAF0(x,y). 

i=l  ^  ^ 

We  can  think  of  this  as  a  surface- fitting  problem  where  the 
data  are  given  by  F^  =  A^F,  i  =  1,2^...,N.  Given  the  data, 
we  can  write  the  approximation  down  iaimedlately . 

We  also  observe  that  if  the  0^  have  support  on  small  sub¬ 
sets  of  D,  and  if  each  also  has  support  on  the  same  set, 

then  the  formula  (6.3)  is  local.  For  example,  if  we  take  A^ 
to  be  point  evaluation  at  the  point  (x^,  y^)  and  0j^(x,  y)  to 
be  a  function  with  support  in  a  neighborhood  of  (x^,yj,),  then 
the  approximation  formula  simply  becomes 
N 

(6.4)  QF(x,y)  =  Zf0  (x,y), 

i=l 

This  is  very  reminiscent  of  the  Lagrange  form  of  interpolation 
(cf.  (4.3)),  but  unless  the  0^^  are  taken  to  satisfy  (4.2), 

QF  will  not  in  fact  be  an  interpolant.  For  this  reason,  for¬ 
mulae  of  the  form  (6.4)  (or  more  generally  (6.3))  are  sometimes 
referred  to  as  quasl-interpolants .  Local  quasi-lnterpolants 
of  the  form  (6.3)  can  be  constructed  simply  by  defining  the 
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N 

functions  with  local  supports.  If  each  of  these  is 

continuous  (or  smooth),  then  QF  will  also  be, 

Although  a  host  of  quasi-interpolants  can  be  constructed 
as  outlined  above,  considerable  care  must  be  exercised  in  order 
to  get  methods  which  give  good  accuracy  (when  the  original 
function  F  is  smooth) ,  As  observed  earlier,  this  is  directly 
related  to  making  the  method  exact  for  polynomials,  i.e.,  such 
that  QP  =  P  for  all  P  in  some  class  of  polynomials. 

To  construct  methods  of  the  form  (6.3)  which  apply  to 

N 

scattered  data,  it  is  necessary  to  construct  appropriate 
While  a  host  of  methods  can  be  constructed  this  way,  it  is  not 
so  easy  to  choose  the  to  make  the  method  exact  for  poly¬ 

nomials  (which,  as  we  remarked  earlier,  is  directly  related  to 
how  well  the  method  will  approximate  smooth  functions  F) .  To 
get  methods  which  do  have  a  reasonable  degree  of  exactness  (and 
a  correspondingly  good  error  bound  for  smooth  functions),  it  is 
easier  to  first  choose  the  {0.}^,  and  then  try  to  find  suit- 
able  While  this  generally  rules  out  using  point  evalu¬ 

ators  at  scattered  data,  it  is  possible  to  construct  methods 
based  on  point  evaluators  at  appropriate  points,  and  such  meth¬ 
ods  can  be  useful  as  second- stage  approximations. 

To  illustrate  these  ideas,  we  consider  construction  of 
local  spline  approximation  methods  following  the  general  treat¬ 
ment  in  Lyche  and  Schuraaker  [131].  Suppose  D  is  enclosed  in 

a  rectangle  H,  and  that  H  is  partitioned  into  subrectangles 

k  Z 

by  a  grid  as  in  (6,2).  Suppose  that  (N. ,  are  the 

Xj  1— X”in^  L*Tl 

tensor  product  B-splines  associated  with  this  partition  (cf. 
(3.30)).  We  are  now  interested  in  approximation  schemes  of  the 
form 

k  £ 

(6.5)  QF(x,y)  =  Z  E  ^'..FN  (x,y). 

i=l-m  j=l-n  J 

In  particular,  we  are  going  to  consider  the  question  of 


constructing  formulae  of  this  type  which  are  exact  for  the  class 

of  polynomials  ,  with  some  fixed  1  <  v  <  m  and  1  <  u  <  n. 

v,u  —  —  —  - 

This  problem  has  a  very  simple  algebraic  solution  if  we  decide 

to  construct  each  A,  .  in  the  form 

ij 


(6.6) 


z  z 

v=l  u=l 


a..  A?.  A^  , 

ijvp  ijv  ijti 


where  the  (A^,  ,  and  (A^,  ,  are  linear  functionals 

ijv  V=1  p=l 

which  apply  to  functions  of  x  and  y  alone,  respectively.  It 
can  be  shown  (cf,  [131])  that  given  any  ijp^ 

tisfying  mild  independence  assumptions,  there  exist  coefficients 

(a.  ,  }  such  that  the  formula  (6.5)  will  be  exact  for  P  . 

ijvp  vu 

In  fact,  these  coefficients  can  easily  be  explicitly  computed. 

To  give  one  example,  suppose 


(6.7) 


(m- 1) 
(n-1) 


i  =  1-m, . . ,,k 

J  —  l”m,  . . . ,  ^ . 


Then  we  obtain 


k  I 

(6.8)  QF(x,y)  =  E  Z  (x,y), 

i=l-m  j=l-n  ^  ^ 

a  formula  which  exactly  reproduces  the  bilinear  polynomials  P. 
This  is  the  multidimensional  (tensor  product)  version  of  the 
Variation  Diminishing  method  of  Marsden  and  Schoenberg;  it  was 
studied  in  some  detail  in  Munteanu  and  Schumaker  [145].  This 
formula  is  closely  related  to  the  Bezier-type  surfaces  construc¬ 
ted  in  Riesenfeld  [163]  (see  also  Gordon  and  Riesenfeld  [89]). 

We  should  observe  that  the  way  formula  (6.5)  now  stands,  it 
may  involve  information  on  F  which  is  taken  from  data  outside 
of  the  domain  D.  This  situation  can  be  rectified  as  follows: 
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Let 

(6.9)  Q  support  D  not  empty}. 

Then  it  can  be  shown  [131]  that  the  method 

(6.10)  QF(x,y)  =  EZ  A  FN  (x,y) 

(i,j)en 


remains  exact  as  long  as  all  functions  are  restricted  to  D. 

To  get  higher-order  methods,  depending  only  on  point  eval¬ 
uations,  we  proceed  as  follows.  Choose 


(6.11) 


X,  <  T.  .  < 

X .  , 

V  =  1,2 

i  ijv 

y.  <  < 

p  =  1,2 

1-m, . . . ,k 

and  j 

=  1-n, . . 

to  be  point  evaluation  at  and  A^ .  to  be  point  evalu- 

^  ijv  IJM 

y 

ation  at  coefficients  in  (6.6)  are  easily  computed. 

Hints  on  where  the  t's  should  be  placed  within  the  support 
of  the  B-splines  are  given  by  the  error  analysis  in  [131]. 

We  close  this  section  with  some  historical  remarks  on  the 


development  of  local  approximation  schemes  in  two  dimensions. 
Early  papers  include  Babuska  [14],  de  Boor  and  Fix  [37],  and 
Fix  and  Strang  [80] .  For  some  methods  involving  triangular 
partitions,  see  Fredrickson  [82].  Quasi-interpolants  were 
constructed  in  de  Boor  and  Fix  [37]  using  point  evaluation 
data,  but  including  derivatives.  We  have  followed  Lyche  and 
Schumaker  [131]  where  general  linear  functionals  are  consider¬ 
ed,  and  where  in  particular,  methods  can  be  constructed  using 
only  point  evaluation  of  the  function.  (Local  integrals  etc. 
would  also  be  possible.)  The  papers  [37]  and  [131]  both  con¬ 
tain  extensive  error  bound  analyses.  It  is  striking  that  these 
local  spline  approximation  methods  give  optimal  order  error 
bounds  for  smooth  functions. 


7 .  Two-stage  processes 


Many  of  the  methods  we  have  discussed  in  this  paper  are 
only  applicable  when  the  data  are  regulany  spaced  (and  in  fact^ 
many  surface- fitting  methods  require  specification  of  derivative 
data  as  well  as  function  values) .  Such  methods  cannot  be  ap¬ 
plied  directly  to  the  scattered  data- fitting  Problem  1,1,  On 
the  other  hand,  some  of  the  most  convenient  local  interpolating 
and  local  approximating  methods  which  do  work  for  scattered 
data  produce  surfaces  which  are  not  globally  smooth  (or  even 
continuous)  ,  Thus,  it  seems  natural  to  consider  the  possibility 
of  constructing  two- stage  processes  in  which  the  first  stage 
uses  the  scattered  data  to  construct  an  approximation  g,  while 
the  second  stage  uses  g  to  generate  data  for  constructing  a 
surface  f  (with  desirable  properties,  such  as  smoothness) , 
Since  it  is  quite  clear  how  various  methods  discussed  in 
the  earlier  sections  might  be  put  together  to  yield  two-stage 
processes,  it  will  suffice  to  mention  just  a  couple  of  examples 
here, 

7,1,  Interpolation/interpolation,  Suppose  that  we  want  to  con¬ 
struct  a  piecewise  bicubic  surface  based  on  data  given  on  a 
rectangular  grid  as  in  (3, 4) -(3, 6).  In  each  subrectangle 
the  16  coefficients  of  the  bicubic  f  (cf,  (4,6))  would  be  de¬ 
termined  by  the  values  of  f,  f  ,  f  ,  and  f  at  each  of  the 

^  X  y  xy 

four  comers.  Now  since  our  original  data- fitting  problem  only 

specifies  the  values  of  the  function  at  the  grid  points,  local 

interpolation  cannot  be  carried  out  directly.  However,  we  can 

use  the  data  to  provide  estimates  for  the  values  of  f  ,  f  , 

X  y 

and  f  at  the  grid  points  (i,e,,  we  construct  g  interpolating 
xy 

the  data);  then  we  can  use  local  bicubic  interpolation  as  a 
second  stage.  The  reader  will  have  no  difficulty  in  imagining 
ways  to  produce  estimates  for  these  quantities.  For  some  meth¬ 
ods  which  appear  in  the  literature,  see  the  papers  of  Akima 
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[5,6],  Koelling  and  Whitten  [121],  and  Spath  [183].  I 

7.2.  Approximation/ interpolation.  Instead  of  making  the  first- 
stage  process  interpolation  as  in  section  7.1,  it  would  also  be 
possible  to  use  an  approximating  process.  For  example,  one 
might  use  least-squares  polynomial  approximation  to  construct 

a  patch  surface  and  then  use  some  convenient  interpolation  pro¬ 
cess  as  a  second  stage.  For  an  example  of  this  type,  see  Hess- 
ing  et  al  [114] . 

7.3.  Approximation/approximation.  This  combination  is  parti¬ 
cularly  convenient  if  we  are  not  concerned  about  getting  an  in¬ 
terpolating  function.  Both  stages  can  be  made  local.  To  give 
an  example,  recently  I  have  constructed  an  algorithm  for  fitting 
surfaces  to  scattered  data  in  which  the  first  stage  consists 

of  polynomial  least-squares  patch  approximation  (with  adaptive 
choice  of  data-- see  section  6),  and  where  the  second  stage  con¬ 
sists  of  direct  local  tensor  product  spline  approximation.  Both 
stages  are  local,  and  the  final  surface  is  a  tensor  product 
spline.  Since  the  second  stage  is  a  direct  method,  it  is  very 
cheap  to  apply.  Experiments  with  real-life  data  (e.g.  from 
heart  potentials,  potential  fields,  and  geological  maps--see 
section  2)  have  produced  very  promising  results.  Details,  in¬ 
cluding  an  analysis  of  error  bounds,  will  appear  elsewhere.  I 
have  also  tried  alternate  versions  where  the  patches  are  con¬ 
structed  as  low-  order  polynomials  which  are  best  approximations 
in  the  or  £  sense  (via  linear  programming)  again  with  adap- 
tive  choice  of  data.  The  results  were  very  similar.  Finally, 

I  have  also  experimented  with  computing  patch  approximations, 
followed  by  continuous  least-squares  tensor-product  spline  ap¬ 
proximation.  Again,  the  experiments  were  promising. 


f 


8.  Contouring 

As  indicated  in  the  introduction,  frequently  the  goal  in 


fitting  a  surface  f  to  data  is  to  construct  a  contour  map  j 
which  approximates  the  contour  map  of  the  unknown  surface  F 
which  produced  the  data.  In  this  section  we  discuss  some  methods 
for  constructing  contour  maps  of  a  surface  f, 

8.1,  Piecewise  linear  functions  on  triangles.  When  the  func¬ 
tion  f  to  be  contoured  is  a  piecewise  linear  function  defined 
on  triangles  (and  globally  continuous),  locating  contours  re¬ 
duces  essentially  to  a  matter  of  good  bookkeeping.  Indeed,  if 
H  is  the  height  of  the  contour  of  interest,  then  it  is  easily 
seen  that  for  a  given  triangle  T  with  vertices.  A,  B,  and  C, 

(8.1)  the  contour  does  not  pass  through  T  if  H  <  min(f(A), 

£(B),  f(C))  or  if  H  >  max(f(A),f(B),f(C)) 

and 

(8.2)  the  contour  intersects  exactly  two  sides  of  T  otherwise. 

If  case  (8.2)  holds,  it  is  easy  to  determine  which  two  sides 
are  intersected  and,  moreover,  by  using  inverse  linear  interpo¬ 
lation  between  vertex  values,  the  points  on  these  sides  where 
the  contour  crosses  can  be  determined.  Specifically,  if,  for 
ex£imple, 

f  (A)  <  H  <  f(B), 

then  the  contour  crosses  the  line  from  A  to  B  at  the  point 
on  the  line  which  is  a  distance  of 


(H-f(A)) 

(f(B)-f(A)) 


B- A 


from  A.  Given  the  points  on  two  sides  of  a  triangle  where  the 
contour  line  crosses,  we  can  now  draw  the  contour  line  since  it 
is  simply  a  straight  line  between  the  points.  An  algorithm  to 
carry  out  this  procedure  requires  enumerating  the  triangles  and 
vertices  and  some  kind  of  effective  search  procedure.  There 
are  several  available  in  the  literature.  For  ALGOL  programs. 
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see  Heap  [111,112],  (An  earlier  paper  of  Heap  and  Pink  [113]  | 

contains  a  similar  FORTRAN  program  but  only  for  regular  triangu- 
larizations .)  Lawson  [128]  discusses  a  similar  algorithm.  The 
algorithms  mentioned  include  two  possible  approaches:  (1)  one 
may  start  with  a  triangle  where  it  is  known  the  contour  inter¬ 
sects,  and  trace  this  contour  as  far  as  it  goes,  or  (2)  one  may 

simply  draw  the  contour  lines  in  all  triangles  which  have  them. 

8.2.  Piecewise  bilinear  functions  on  rectangles.  Suppose  now 
that  the  function  f  to  be  contoured  is  a  piecewise  (continuous) 
function  on  a  rectangle  partitioned  into  subrectangles  by  a  grid. 
Since  f  is  linear  in  x  or  y  on  the  edges,  it  follows  that 

we  can  again  determine  whether  a  contour  line  of  height  H 

crosses  an  edge  by  inverse  linear  interpolation.  There  is  in 
this  case,  however,  a  serious  difficulty  which  does  not  arise 
in  the  case  of  triangles.  It  may  happen  that  the  height  H 
lies  on  three  or  even  four  sides  of  the  rectangle.  In  this 
case,  it  is  possible  that  two  different  contour  lines  pass 
through  the  rectangle,  and  it  is  not  clear  how  to  interconnect 
the  points  (see  Figure  6) . 


Figure  6,  Two  Contours  in  a  Rectangle 
Put  another  way,  if  we  are  following  a  contour  and  enter  a  rec¬ 
tangle  as  shown  above  in  Figure  6  on  the  bottom  line,  then  it 
is  not  clear  whether  we  should  now  turn  right  or  turn  left.  One 
approach  to  designing  an  algorithm  in  this  case  is  to  simply 

I 

( 


always  go  right,  say,  even  though  this  may  in  the  end  be  wrong. 
(If  it  is,  we  have  to  start  over  with  a  coarser  mesh.)  This 
technique  was  incorporated  in  an  algorithm  by  Heap  [111,112]-- 
the  paper  contains  a  FORTRAN  program,  (An  earlier  ALGOL  pro¬ 
gram  can  be  found  in  Heap  and  Pink  [113]. 

A  second  approach  to  handling  the  ambiguity  problem  is  com 
pute  an  approximation  to  the  value  of  f  at  the  center  of  the 
rectangle  (e.g.,  by  taking  the  average  of  the  four-comer  val¬ 
ues)  and  then  to  triangulate  the  rectangle.  This  amounts  to  a 
second-stage  approximation  process,  and  the  surface  contoured 
is  no  longer  f  itself  but  an  approximation  g.  This  idea  was 
programmed  in  ALGOL  in  Heap  and  Pink  [113]  and  in  FORTRAN  in 
Heap  [111,112]. 

Once  the  set  of  points  for  a  particular  contour  have  been 
found,  there  are  a  variety  of  ways  of  drawing  a  contour  line 
through  these  points.  One  possibility  is  to  simply  draw 
straight  lines  between  each  of  the  points.  The  actual  contour 
lines  are  expressions  of  the  form  y  =  (a+bx) / (c+dx)  in  each 
rectangle.  These  are  generally  not  straight  lines.  Hence,  if 
smoother  contours  are  desired,  one.  may  use  any  one  of  a  number 
of  methods  for  drawing  a  smooth  curve  through  an  ordered  set  of 
points  in  the  plane.  For  example,  the  curve  could  be  computed 
in  parametric  form  using  one- dimensional  splines.  Another  pos¬ 
sibility  would  be  to  use  the  Bezier  methods  with  either  Bern¬ 
stein  polynomials  or  with  B-splines  (cf.  Gordon  and  Riesenfeld 
[89]  and  Riesenfeld  [163]),  although  in  this  case  the  curves 
will  not  exactly  go  through  the  points.  For  other  algorithms 
see  Marlow  and  Powell  [138]  or  McConalogue  [139]. 

8.3.  Piecewise  quadratics  on  triangles.  Suppose  now  that  f 
is  a  piecewise  quadratic  defined  on  a  triangular  partition.  In 
this  case  a  contour  line  at  height  H  passing  through  a  trian¬ 
gle  must  be  described  by  a  conic  section.  Such  a  section  can 


4 


49 


be  represented  in  parametric  form  as 

X(t)  =  (bQ4b^t+b2t^)/(b^4b^t  )  b^tS 

y(t)  =  (b^ +b^t  ^  bgt^)/(b^  ^  b^t  +  b^t^), 

see  Powell  [156].  Powell  has  promised  an  algorithm  based  on 
this  observation. 

We  turn  now  to  some  methods  for  handling  general  functions 
f  on  arbitrary  domains  D. 


8.4.  A  simple  line-printer  method.  The  following  simple-minded 
method  can  produce  reasonable- looking  contours  without  excessive 
computation,  and  without  recourse  to  a  p’.occor.  Suppose  H  is 
a  rectangle  enclosing  the  domain  D,  and  that  we  partition  H 
as  H  =  with  a  rectangular  grid  as  in  (6.2).  Let  HL  <  HU 

be  given  real  numbers.  Finally,  suppose  that  t^^  is  some 
point  in  where  f  can  be  evaluated  (perhaps  one  of  the 

comers  or  the  center) .  Define 

r 


(8.3) 


'ij 


0 

=  <  9 


if  f(t_,  .)  <  HL 
ij 


if 


f(t^.)  >HU 


V  ,  if  HL  +  (v-l)h  <  f(t^j)  <  HL  +  vh,  l<v<8. 


for  all  i  =  0, 1, . .  .,k  and  j  =0,1,...,!  (where  h  =  (HU-HL)  /8) . 
The  (k+2)  by  (i+2)  matrix  C  contains  only  integers,  and  if  it 
is  printed  out  without  either  horizontal  or  vertical  spacing, 
we  obtain  a  reasonable- looking  contour  map  of  the  function.  A 
typical  example  is  included  in  Figure  7.  The  method  can  be 
refined  by  using  an  alpha-numeric  array  C  and  more  than  10 
symbols.  It  can  also  be  refined  by  using  a  printer  with  appro¬ 
priate  horizontal  spacing  so  that  each  symbol  occupies  a  square 
rather  than  a  rectangle  (e.g.,  cf.  Buneman  [40]). 

8.5.  Threading  on  a  rectangular  grid.  As  in  section  8.4, 
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Figure  7.  A  Simple  Contour  Map  (Heart  Potential) 


I  suppose  that  D  is  imbedded  in  a  rectangle  H  which  has  been  | 
partitioned  by  a  rectangular  grid  as  in  (6.2).  Assuming  that 
f  is  continuous^  it  is  still  possible  to  decide  which  of  the 
grid  lines  a  particular  contour  of  height  H  crosses  by  examin¬ 
ing  the  end-points  of  each  such  line.  Since  f  is  not  generally 
linear  along  such  a  line,  we  cannot  determine  exactly  where  the 
crossing  point  is  by  linear  inverse  interpolation.  However,  if 
we  are  willing  to  evaluate  f  a  few  times  along  this  line,  we 
can  estimate  the  crossing  point  quite  accurately  by  bisection. 


for  example.  Once  a  sequence  of  points  on  a  contour  has  been 
determined,  we  may  thread  a  curve  through  the  points  just  as  in 
section  8.2. 

This  method  does  have  one  serious  drawback,  however, -- 
just  as  with  the  method  discussed  in  section  8.2--,  if  we  are 
tracing  a  contour  it  may  happen  that  after  entering  a  triangle 
there  is  an  ambiguity  as  to  which  of  two  points  to  use  to  exit 
the  rectangle.  One  could  opt  for  an  ad  hoc  rule  or  try  the 
second- stage  approximation  described  in  section  8,2.  For  an 
example  of  how  this  method  works,  see  Falconer  [76]  (based  on 
Lodwick  and  Whittle  [130]),  where  it  is  applied  to  a  surface 
constructed  by  local  weighted  quadratic  polynomial  least-squares 
approximation.  Since  bisection  is  involved,  one  should  realize 
that  in  drawing  contours  with  this  routing  the  surface  f  is 
going  to  be  evaluated  a  great  many  times. 

8.6.  Threading  on  a  triangular  grid.  An  obvious  cure  for  the 
ambiguity  discussed  in  section  8,5  for  threading  on  a  rectangu¬ 
lar  grid  is  to  use  a  triangular  partition  in  the  first  place. 
Then  the  bisection  method  coupled  with  a  threading  routine  leads 
immediately  to  a  contouring  routine  for  general  surfaces  f. 
Strangely  enough,  I  have  not  been  able  to  find  anywhere  where 
this  method  has  been  suggested. 

I  have  made  no  effort  to  track  down  all  available  papers 
on  contouring.  A  few  which  I  did  find  and  have  not  yet  men¬ 
tioned  are  Cottafawa  and  le  Moli  [60],  Dayhoff  [64],  and  Pelto 
et  al  [152].  There  are  many  others. 

In  some  cases  it  may  be  desirable  to  have  a  more  graphic 
picture  of  a  surface  than  a  contour  map  can  provide.  Recently 
there  has  been  considerable  effort  devoted  to  computer  methods 
for  displaying  surfaces  on  a  scope  or  with  a  plotter.  For  some 
examples  of  output  and  a  discussion  of  methods,  see  e.g.  the 
book  by  Barnhill  and  Riesenfeld  [20]  on  computer-aided  design. 


If  an  actual  3-D  picture  is  desired  instead  of  just  a  perspec-  | 
tive,  it  is  even  possible  to  produce  holographs. 
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